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1.  INTRODUCTION 


Our  research  objectives  are  directed  toward  understanding  electromagnetic  propagation  in 
dielectric  media  with  conductive  losses;  investigation  of  the  complementary  use  of  both 
seismic  and  ground  penetrating  radar  (GPR)  technologies;  investigating  the  use  of  lattice- 
gas  cellular  automata  methods  to  determine  the  characteristics  of  the  ground  water 
hydrology  from  non-invasive  measurements  made  at  the  ground  surface;  the  development 
of  algorithms  and  approaches  to  create  a  new  improved  5x5  global  grid  of  Rayleigh  and 
Love  surface  wave  group  velocities;  and  the  collection  and  interpretation  of  seismic  data 
to  determine  crustal  structure.  Oftentimes  to  properly  evaluate  a  research  application,  it 
has  to  be  tested  with  new  and  old  seismic  data  bases  and  these  results  compared  to  those 
using  existing  approaches.  Papers  and  abstracts  of  papers  giving  more  detail  on  these 
efforts  are  included  in  the  Appendix. 

2.  ELECTROMAGNETIC  PROPAGATION  IN  DIELECTRIC  MEDIA 
WITH  CONDUCTIVE  LOSSES 

Understanding  the  propagation  of  radar  signals  in  dielectric  media  with  conductive  losses 
is  fundamental  to  GPR  measurement  interpretation.  With  a  matrix  method  which 
describes  plane-wave  propagation,  we  have  attempted  to  model  electromagnetic  and 
elastic  waves  in  layered  media.  The  Ray  theory  (Goodman,  1994)  and  finite  difference 
techniques  (Roberts,  1994)  have  been  used  to  model  electromagnetic  waves  in  the  earth; 
however  dispersion  and  attenuation  are  non-trivial  to  incorporate  into  these  time-domain 
algorithms.  The  main  advantages  of  the  matrix  method  are:  (1)  it  is  an  exact  solution; 
(2)  it  is  developed  in  the  frequency  domain  so  that  attenuation  and  dispersion  effects  are 
easily  incorporated;  (3)  results  are  obtained  in  real-time  on  a  workstation,  without  the 
requirement  of  a  supercomputer  for  the  calculations;  and  (4)  extensions  of  the  method 
model  point  scatterers  and  lateral  heterogeneity  in  the  earth  media. 

The  matrix  method  for  solving  electromagnetic  wave  propagation  outlined  above  is  more 
generic  than  simply  a  solution  to  Maxwell’s  equation  (Ursin,  1983).  As  a  solution  to  an 
eigenvalue  problem,  the  matrix  method  has  been  used  to  model  elastic  waves  in  a  realistic 
earth  (Kennett  et  al.,  1979);  e.g.,  we  have  used  it  to  model  compressional  and  shear  waves 
in  anisotropic  media.  The  ability  to  model  both  radar  and  seismic  data  will  greatly 
benefit  the  effort  toward  combining  the  two  technologies. 
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Seismic  waves  are  sensitive  to  the  mechanical  properties  of  the  shallow  subsurface,  and 
radar  is  sensitive  to  the  chemistry  in  the  soils.  In  the  course  of  research,  we  explored  this 
difference  between  the  data  sets,  and  worked  on  understanding  more  of  the 
complementary  nature  of  the  two  non-invasive  geophysical  methods. 

For  synthesis  of  GPR,  we  arranged  Maxwell’s  equations  describing  the  physics  of 
electromagnetic  fields  in  general  media  into  a  system  of  equations  to  be  solved  as  an 
eigenvalue  problem.  A  differential  system  describes  the  propagation  of  the  waves  within 
the  plane-wave  approximation.  The  local  eigenvalues  of  the  system  are  the  permissible 
vertical  wavenumbers  and  include  attenuation,  while  the  local  eigenvectors  are  used  to 
match  boundary  conditions  across  interfaces  between  layers  in  the  medium.  The 
eigenvectors  and  eigenvalues  are  incorporated  into  a  fundamental  matrix  describing 
propagating  waves  in  the  realistic  earth,  where  vertical  variations  in  permittivity, 
permeability  and  conductivity  (as  functions  of  frequency)  can  all  be  modeled 

The  fundamental  matrix  which  describes  the  propagation  of  electromagnetic  energy 
through  a  layered  medium  gives  exactly  the  reflection  and  transmission  response  of  the 
medium  for  a  monochromatic  plane  wave  at  a  given  angle  of  incidence.  Calculation  of 
the  amplitude  variation  has  direct  applicability  to  depth  of  penetration  (skin-depth) 
studies;  however,  there  is  more  potential  information  in  the  calculation  since  the  matrix 
method  is  a  full  wave  optics  solution,  where  amplitude  and  phase  are  calculated  together 
(as  opposed  to  ray  optics  in  which  amplitude  and  phase  are  separate  computational 
problems).  We  have  explored  ways  to  use  the  phase  information  in  the  GPR  signal  for 
detection  of  changes  in  physical  quantities. 

Calculation  of  the  reflection  response  at  many  different  frequencies,  and  integrating  over 
those  frequencies,  gives  synthetic  radar  signals  in  the  time  domain.  We  used  a  narrow- 
band  zero-phase  wavelet  (a  Ricker  wavelet)  as  our  source  function  for  modeling  radar 
data.  The  matrix  method  deals  directly  with  the  vector  nature  of  the  electric  and  magnetic 
fields,  so  that  we  can  explore  source  and  medium  polarization  effects.  The  most  common 
form  of  ground  penetrating  radar  data  collection  is  in  monostatic  mode,  which  uses  the 
geometry  of  a  coincident  source  and  receiver  and  yields  vertical  incidence  radar  returns. 
Calculation  of  the  time-domain  response  at  zero  horizontal  wavenumber  generates  a 
vertical-incidence  synthetic  radar  record  and  for  more  than  one  receiver  gives  a  synthetic 
monostatic  radar  record  section.  We  introduced  lateral  variation  in  the  earth  model  from 
one  receiver  to  the  next  in  order  to  approximate  monostatic  record  sections  over 
structures  that  are  not  horizontally  planar. 
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The  eigenvalue  solution  to  wave  propagation  requires  the  medium  be  horizontally 
homogeneous.  Introducing  a  perturbation  to  the  laterally  homogeneous  background 
model  is  non-trivial,  mainly  in  that  it  requires  a  convolution  involving  the  response  to  the 
perturbation  and  the  full  wavefield  itself.  With  a  point  scatterer  as  the  perturbation  to  the 
background  earth  model,  and  assuming  near- vertical  incidence  waves  upon  the  scatterer, 
we  hope  to  achieve  analytic  results  for  the  convolution  integral,  giving  an  approximate 
form  for  the  3-D  response  of  a  scatterer  that  can  be  introduced  into  the  1-D  matrix 
algorithm.  With  each  scatterer  acting  as  a  'Huygens’  source,  a  superposition  of  the 
response  from  numerous  scattering  points  can  simulate  laterally  heterogeneous  media. 

3.  FIELD  TEST  OF  THE  COMPLEMENTARY  TECHNOLOGIES: 

SEISMIC  AND  GPR 

Coupled  with  the  computer  modeling  effort,  there  were  field  efforts  in  which  we 
coordinated  the  collection  of  both  ground  penetrating  radar  and  shallow  high  resolution 
seismic  data  at  one  or  more  sites.  The  data  collected  at  these  sites  has  allowed  us  to 
evaluate  our  near-surface  Earth  model  and  continue  to  fine  tune  it,  when  necessary,  for 
the  individual  sites. 

Our  researchers  and  programmers  are  involved  with  the  processing,  analyzing  and 
interpretation  of  geophysical  data  from  field  experiments  at  Dover  AFB  in  1995  and  1996 
and  at  Hanscom  AFB  in  1996.  The  1996  field  tests  collected  only  seismic  refraction  data 
as  a  check  on  seismic  velocity  structure  and  to  'image'  the  water  table,  something  neither 
of  the  1995  instruments  was  able  to  measure.  The  1995  Dover  data  base  consists  of 
shallow  seismic  reflection  records,  ground  penetrating  radar,  and  electromagnetic 
induction  measurements.  Additionally,  data  taken  by  other  researchers  at  Dover  AFB  in 
1995  was  exchanged.  We  created  a  data  base  of  all  the  measurements  collected,  and 
distributed  it  to  all  the  researchers.  This  data  base  includes  terrain  conductivity,  elevation 
information,  the  seismic  and  radar  data  we  collected,  and  the  radar  data  collected  by  other 
researchers.  Initial  analysis  of  these  measurements  has  been  completed.  This  analysis 
provided:  (1)  comparison  of  seismic  sources  tested  at  Dover  AFB  for  common  mid-point 
reflection  profiling  and  seismic  characterization  of  the  Groundwater  Remediation  Field 
Lab  (GRFL)  site;  (2)  description  of  the  seismic  data  and  the  low  frequency  ground 
penetrating  radar  acquired  at  Dover  AFB  as  one  of  the  first  direct  comparisons  of 
coincident  shallow  high  resolution  seismic  data  and  ground  penetrating  radar  that  image 
similar  targets;  (3)  comparison  of  ground  penetrating  radar  collected  by  the  two  most 
commonly  used  systems  with  an  array  of  source  frequencies  for  each  instrument;  and  (4) 
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integration  of  all  geophysical  site  characterization  data  at  the  GRFL  as  an  extensive  case 
study  for  a  well  characterized  shallow  aquifer  site. 

Our  scientists  organized  a  simple  field  experiment  at  Hanscom  AFB  in  order  to  test  field 
equipment  that  was  to  be  used  in  a  field  experiment  at  Dover  AFB.  As  part  of  the 
Hanscom  experiment,  they  conducted  a  seismic  refraction  survey  near  to  a  contaminated 
site  on  Hanscom  Air  Field.  Initial  data  analysis  was  meant  to  help  determine  field 
geometry  and  data  acquisition  parameters  for  another  field  experiment  at  Dover  AFB. 
Further  work  will  involve  determination  of  bedrock  profile  map  under  seismic  acquisition 
line  at  the  Hanscom  site.  Our  researchers  organized  and  help  conduct  a  field  experiment 
at  Dover  AFB,  Dover  DE.  We  conducted  a  seismic  refraction  survey  of  the  Groundwater 
Remediation  Field  Laboratory  site.  The  primary  goals  were  to  map  —if  possible-  the 
aquitard  boundary  (the  confining  clay  unit  for  the  perched  aquifer  at  the  site)  and 
determine  general  seismic  velocity  structure  of  the  shallow  subsurface  to  complement 
seismic  reflection  work  carried  out  in  1995.  Initial  analysis  of  the  data  has  been  directed 
toward  the  second  of  these  two  goals. 

4.  USE  OF  LATTICE-GAS  CELLULAR  AUTOMATA  METHODS 

Cellular  automata  form  a  subclass  of  artificial  neural  networks  which  can  be  described  as 
lattices  of  individual  finite  state  automata  updated  in  discrete  time  steps,  with  uniform 
structure  and  polynomial  connectivity.  The  state  of  each  lattice  site  at  a  given  time  step  is 
defined  in  a  deterministic  or  probabilistic  way  by  a  given  neighborhood  of  lattice  sites  at 
a  previous  time  step. 

Special  cases  of  cellular  automata  are  capable  of  describing  complex  collective  behavior. 
The  simulation  of  fluid  dynamics  is  an  example  which  proves  cellular  automata  are 
complementary  and  powerful  tools  to  model  phenomena  that  would  normally  be  the 
exclusive  domain  of  partial  differential  equations. 

This  class  was  named  lattice-gas  cellular  automata  and  has  not  only  become  a  ’toy  model' 
for  the  exploration  of  the  microscopic  basis  of  hydrodynamics,  but  also  tools  for  the 
numerical  study  of  complex  problems  in  fluid  mechanics. 

While  the  lattice  gas  may,  in  principle,  be  used  for  nearly  any  problem  in  hydrodynamic 
simulation,  many  of  the  most  successful  applications  have  involved  either  a  complex 
fluid,  a  complex  geometry,  or  both. 
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Lattice-gas  cellular  automata  methods  (LGCA)  are  an  attempt  to  use  the  theory  of 
discrete  molecular  dynamics  to  model  physical  processes  such  as  fluid  flow  through 
porous  media  and  wave  propagation  in  heterogeneous  media.  Solving  these  problems  via 
separate  partial  differential  equations  keeps  the  parameterizations  separate.  Using  the 
single  LGCA  theory  to  describe  both  of  these  parameters  (e.g.,  elastic  wave  velocities)  to 
the  desired  physical  properties  (e.g.,  hydraulic  conductivity).  This  type  of  mapping  is 
most  important  in  the  environmental  area  where  it  is  often  necessary  to  determine  the 
characteristics  of  the  ground  water  hydrology  from  non-invasive  geophysical 
measurements  made  at  the  ground  surface. 

5.  SURFACE  WAVE  REGIONALIZATION 

Our  scientists  and  researchers  have  attempted  to  improve  the  existing  global  grids  of 
Rayleigh  and  Love  surface  wave  group  velocities.  The  grid  is  used  to  estimate  surface 
wave  arrival  times  for  use  in  seismic  surface  wave  association  codes;  in  other  words,  for 
determining  which  of  many  recorded  surface  waves  is  associated  with  a  particular  seismic 
event.  In  the  original  grid,  the  emphasis  was  on  surface  wave  energy  in  the  17-23  second 
band  used  for  measuring  the  surface  wave  magnitude,  Ms  which  is  used  for  yield 
estimation  and  as  a  discriminant  when  compared  with  the  body  wave  magnitude,  mb- 
Additionally,  this  grid  contained  16  Rayleigh  wave  regions  of  which  9  were  oceans  or 
shallow  seas;  South  America,  Australia,  and  Antarctica  were  considered  the  same  group 
velocity  region.  The  drawback  with  this  model,  excluding  the  17-23  second  range,  the 
velocities  and  more  important  the  shape  of  these  regionalized  dispersion  curves  from  10- 
50  seconds  did  not  agree  with  observations  and  were  of  no  use  for  phase  matched 
filtering,  a  common  technique  for  finding  and  extracting  surface  wave  signals  from  noise. 

As  a  starting  point,  we  requested  and  received  from  the  Air  Force  Technical  Applications 
Center  (AFT AC)  their  regional  surface  wave  dispersion  subroutines  and  data  base,  so  that 
we  could  include  new  global  crustal  structure  information  in  it.  The  global  model  we 
initially  used  was  Mooney's  1994  Global  Crustal  Model.  This  model  is  described  on  a 
5x5  degree  grid  and  gives  the  thickness  of  ice,  water,  soft  sediments,  hard  sediments, 
upper  crust,  and  lower  crust.  The  initial  version  we  received  additionally  had  an  eight 
layer  P-velocity  model  of  the  regions.  We  added  shear  velocity  and  density  to  Mooney's 
model.  However,  we  later  received  Mooney,  Laske  and  Masters  (1995)  Global  Crustal 
Models  with  gridded  global  crustal  structure.  This  model  was  an  improvement  over  the 
1994  model  in  that  there  are  now  89  distinct  structure  types  versus  the  previous  52  types 
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and  shear  velocities  and  densities  have  added  the  8  layers  making  up  each  of  the  89 
models.  Our  researchers  then  incorporated  the  models  into  AFT AC's  codes  for  surface 
wave  ray  or  path  tracing.  We  completed  the  calculation  of  the  fundamental  Rayleigh  and 
Love  phase  and  group  velocities  for  the  1995  Global  Crustal  Models  including  the  crust 
mantle  interface  (Moho)  at  a  5x5  degree  grid  or  each  crust  type  at  specified  periods  using 
our  dispersion  code.  The  periods  of  AFTAC's  observed  velocities  were  overlapped  by 
starting  at  38  seconds  and  going  down  to  7  seconds.  Global  contour  maps  of  Rayleigh 
and  Love  group  velocities  at  periods  of  20  and  10  seconds  were  plotted  Our  original 
purpose  for  using  crustal  structures  was  for  calculating  regional  dispersion,  i.e.,  for 
periods  less  than  20  seconds. 

Using  over  4000  source-receiver  pairs,  Stevens  (1996)  was  able  to  perform  a  tomographic 
determination  of  Rayleigh  group  velocities  from  300  to  12.5  seconds  and  refine  the  Okal 
dispersion  values  in  the  shield  and  sub-divide  the  mountain  region  into  two.  In  addition, 
he  obtained  vastly  improved  group  velocities  for  the  trench  and  subduction  region  for  a 
total  of  8  distinct  dispersion  regions.  This  10x10  degree  model  is  implemented  at  the 
Prototype  International  Data  Center.  Again  using  Global  Tectonic  Regionalization 
(GTR1)  as  a  guide,  we  interpolated  the  Stevens-Okal  regionalization  into  a  5x5  degree 
global  model.  Arbitrarily  reducing  the  grid  size  for  tomographic  determined  velocities 
can  create  problems,  but  in  this  case  the  inversion  path  lengths  across  the  age  boundaries 
for  the  oceanic  determined  dispersion  curves  of  Yu  &  Mitchell  (1979)  and  Mitchell  &  Yu 
(1980)  were  based  on  smooth  boundaries,  and  so  in  the  oceans  the  smaller  the  grid  the 
closer  to  the  actual  path  lengths  used  in  the  original  inversions.  This  was  also  added  to 
our  composite  model. 

We  received  test  data  from  AFTAC  which  was  to  be  used  compare  our  model’s  to  five 
other  current  models.  The  test  data  was  for  five  events;  three  events  are  in  the  Nevada 
Test  Site  (NTS),  one  in  China,  and  one  in  Kazakh.  Each  data  set  had  the  location  of  the 
stations  recording  the  events  and  the  observed  travel  time  and  great  circle  group 
velocities  to  the  stations.  We  were  able  to  read  the  test  data  and  calculate  a  group  travel 
time,  between  the  events  and  recording  stations  in  the  test  data  file  for  each  of  the  six 
proposed  regionalizations  and  their  regional  group  velocity  dispersions  curves.  Statistics 
were  determined  for  all  five  events  such  as  the  number  of  stations,  the  distribution  of 
range  and  azimuthal,  and  the  average  period  for  which  a  group  arrival  time  was 
determined.  For  the  five  events  tested,  the  composite  model  was  the  best  in  reduction  of 
the  average  difference  over  the  stations  and  it's  standard  deviation. 
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6.  COLLECTION  AND  INTERPRETATION  OF  SEISMIC  DATA 
TO  DETERMINE  CRUSTAL  STRUCTURE 

In  an  effort  to  determine  the  crustal  and  upper  mantle  structure  in  a  wide  region  around 
the  Pinedale  Seismic  Array.  Our  researchers  helped  install  a  large-aperture  seismic  array 
surrounding  the  existing  Pinedale  Seismic  Array  to  collect  data  for  the  purpose  of 
characterizing  the  region  geophysically.  The  array  operated  for  2.5  months  starting  near 
the  end  of  July  1995  and  consisted  of  five  broadband  stations  each  located  approximately 
80-100  km  from  the  Pinedale  Seismic  Array.  Each  station  contains  broadband  and/or 
long-period  sensors  installed  in  a  shelter  mounted  to  a  concrete  pad,  a  refraction 
technology  data  acquisition  system,  a  530-Mbyte  hard  disk,  and  a  GPS  clock  for  accurate 
timing.  Power  at  each  site  came  from  solar  panels  that  we  installed  or  by  direct  line 
power.  The  large-aperture  array  was  designed  to  provide  data  for  three  separate  studies: 

•  measure  teleseismic  surface  wave  propagation  (phase  and  group 
velocity,  attenuation)  across  the  region, 

•  receiver  function  analysis  of  teleseismic  P-waves  at  individual  array  stations, 

•  modeling  of  surface  and  body  waves  from  regional  earthquakes. 

Additionally,  we  recorded  and  are  interpreting  a  pair  of  large  chemical  explosions 
detonated  in  the  western  United  States  in  August-September  1995.  These  two 
Ammonium  Nitrate  Fuel  Oil  (ANFO)  explosions  were  detonated  within  a  week  of  each 
other  near  the  Gas  Hills,  Wyoming  for  the  purpose  of  investigating  the  structure  of  the 
crust  and  upper  mantle.  These  explosions  were  monitored  at  varying  seismic  sites  along 
the  Green  River  Seismic  Profile,  including  one  site  that  was  common  to  both  explosions. 
The  data  were  collected  and  processed  at  the  common  site  for  both  explosions.  These 
explosions  were  recorded  by  the  Pinedale  Large-Aperture  Array  installed  around  the 
existing  seismic  array.  The  second  source  of  measurements,  used  to  described  the 
refraction  profile,  came  from  47  individual  seismic  stations  that  were  installed  and 
deployed  for  a  period  up  to  three  days.  The  stations  were  uniformly  spaced  by  3.2  km 
starting  from  Big  Sandy,  Wyoming,  west  to  the  border  of  Wyoming  and  Idaho.  The  data 
from  the  47  stations  were  displayed  for  two  chemical  explosions  detonated  near  the  Gas 
Hills,  Wyoming.  The  profile  data  collected  is  essential  to  determine  the  crustal  structure 
beneath  the  Wind  River  Mountains.  Included  in  the  47  stations  were  three  different  types 
of  station  configurations.  The  data  was  recorded  on  three  different  instruments.  These 
recorders  are  Terra  Technology  Recorders,  Refraction  Technology  Recorders,  and  PDAS 
Recorders,  each  having  a  different  format  for  storing  data.  Therefore,  the  data  had  to  be 
processed  in  a  different  manner  to  produce  a  uniform  data  set.  This  involves  numerous 
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quality  checks  of  the  data  before  and  after  processing.  The  processing  of  the  seismic  data 
in  a  similar  fashion  as  the  other  station  configurations  allowed  us  to  combine  the  data  and 
display  it  as  a  single  refraction  profile. 

The  refraction  data  set,  including  the  two  blasts  recorded  at  the  47  seismic  stations,  was 
downloaded  onto  an  Exabyte  tape  which  is  now  available  to  other  scientists  who  want  to 
use  the  Wyoming  refraction  data  set.  More  analysis  was  conducted  on  the  similarity  of 
the  two  ANFO  explosions  used  in  this  experiment.  The  data  for  both  explosions  were 
collected  and  processed  at  a  common  site.  Spectral  amplitude  ratios  were  examined  for 
different  seismic  phases  including  a  pre-event  noise  sample. 
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PRELIMINARY  GEOPHYSICAL  CHARACTERIZATION  OF  GROUNDWATER 
REMEDIATION  FIELD  LABORATORY,  DOVER  AIR  FORCE  BASE,  DELAWARE 


Kadinsky-Cade,  Katherine 

Phillips  Laboratory  Geophysics  Directorate,  Earth  Sciences  Division 

Cardimona,  Steven 

Institute  for  Scientific  Research,  Boston  College 


ABSTRACT: 

During  a  two  week  period  in  June  1995,  the  Earth  Sciences  Division  of  Phillips  Laboratory 
performed  an  integrated  geophysical  survey  at  a  site  that  will  be  used  for  a  series  of 
groundwater  remeSiationstudies.  This  survey  was  carried  out  with  neld  support  by  iiic 
Colorado  School  of  Mines,  Kansas  Geological  Survey  and  Elohi  Geophysics,  Inc.  The 
groundwater  remediation  work  will  be  part  of  the  Air  Force’s  contribution  to  the  Strategic 
Environmental  Research  and  Development  Program  (SERDP).  The  geophysical  survey 
helps  provide  site  characterization  information  that  can  aid  in  the  siting  of  test  cells  for 
future  remediation  experiments.  An  important  focus  of  this  work  is  to  compare  coincident 
surface  seismic  and  ground  penetrating  radar  (GPR)  measurements.  The  site  consists  of  a 
10  to  15  meter  thick  sandy  aquifer  overlying  a  clay  aquitard.  The  site  is  heterogeneous, 
with  clay  lenses  and  gravel  present  within  the  aquifer.  A  variety  of  seismic  sources  were 
tested  at  the  site,  including  a  Betsy  firing  rod  shooting  12  gauge  blanks,  a  sledgehammer,  a 
slidehammer  and  a  portable  vibrator.  We  used  100  Hz  geophones  throughout  the  seismic 
survey.  Recording  was  done  in  a  common  depth  point  reflection  profiling  configuration. 
A  downhole  seismic  survey  was  performed  by  lowering  a  hydrophone  into  a  monitoring 
well  at  the  site.  Continuously  recorded  GPR  profiles  at  300  and  500  MHz  were  collected 
to  complement  a  set  of  100  and  200  MHz  measurements  obtained  by  the  University  of 
Delaware  at  the  same  site.  The  surface  data  are  being  supplemental  by  cone 
penetrometer  measurements  (resistivity  and  uphole  seismic)  performed  by  Applied 
Research  Associates,  Inc.  We  made  terrain  conductivity  measurements  to  aid  in  the 
interpretation  of  the  GPR  data.  Precise  terrain  elevation  measurements  are  being  used  for 
static  corrections.  Data  processing  is  currently  in  progress.  It  is  commonly  pointed  out 
that  the  seismic  and  GPR  techniques  tend  to  be  successful  in  mutually  exclusive  areas  (GPR 
attenuated  by  high  conductivity  clay  media  and  seismic  plagued  by  poor  coupling  in 
sandy  soils).  At  this  site  both  techniques  seem  to  work  reasonably  well,  and  variations  in 
the  imaging  capability  of  each  technique  across  the  site  will  be  the  subject  of  an  interesting 
research  study. 
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MODELING  GROUND  PENETRATING  RADAR  TRANSMISSION  PROPERTIES  IN 
CONDUCTIVE  AND  DISPERSIVE  MEDIA 


Steve  Cardimona.  Institute  for  Scientific  Research,  Boston  College 
Katharine  Kadinsky-Cade,  Air  Force  Phillips  Lab,  Earth  Sciences 
Both  at:  PL/GPE,  29  Randolph  Road,  Hanscom  AFB,  MA,  01731-3010 

Maxwell’s  equations  describing  the  propagation  of  electromagnetic  plane  waves  in  dielectric  media  with  conductive 
losses  are  solved  exactly  by  finding  solutions  to  an  eigenvalue  problem  for  electric  and  magnetic  field  intensities  in 
media  with  depth-dependent  properties.  For  a  given  frequency  and  horizontal  wavenumber  (propagation  angle),  the 
local  eigenvalues  of  the  differential  system  define  the  permissible  vertical  wavenumbers  within  the  medium  and 
include  the  attenuation  due  to  non-zero  conductivity.  The  matrix  of  eigenvectors  and  its  inverse  are  used  to  match 
the  boundary  conditions,  i.e.,  the  continuity  of  the  tangential  components  of  the  electric  and  magnetic  field 
intensities  across  the  interfaces  between  layers.  Using  the  local  eigenvalues  and  eigenvectors  of  the  system  we 
compute  a  fundamental  matrix  for  the  layered  medium,  from  which  we  obtain  the  reflection  and  transmission 
properties  of  the  complete  subsurface  structure.  The  time-domain  wavefield  solution  can  then  be  built  up  through  a 
summation  (slow-Fourier  transform)  over  the  contributions  at  all  horizontal  wavenumbers  of  interest  and  an 
integration  (a  fast  Fourier  transform)  over  the  frequencies  in  the  source  wavelet.  The  reflection  coefficient  can  be 
integrated,  for  example,  to  yield  the  backscattered  wavefield  which  allows  us  to  synthesize  surface  collected  ground 
penetrating  radar  data.  We  use  this  method  to  explore  the  transmission  properties  of  attenuative  and  dispersive 
layered  media  by  varying  the  parameters  of  conductivity,  magnetic  permeability,  electric  permittivity,  peak  probing 
frequency,  and  angle  of  incidence.  Computation  of  the  transmissivity  as  a  function  of  depth  for  a  variety  of 
parameters  is  accomplished  rapidly  since  the  integration  over  wavenumber  and  frequency  to  get  the  time-domain 
wavefield  is  not  performed.  In  this  way  we  can  easily  investigate  the  sensitivity  of  the  ground  penetrating  radar 
signal  to  variations  in  realistic,  near-surface  medium  properties. 
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Isotropic  and  Deviatoric  Moment  Inversion  of  Regional  Surface  Waves 
from  Nevada  Test  Site  Explosions: 

Implications  for  Yield  Estimation  and  Seismic  Discrimination 

Bradley  B.  Woods  and  David.  G.  Harkrider 
Seismological  Laboratory,  CalTech,  Pasadena,  CA 
AFOSR  F49620-93- 1-0221 

Institute  of  Scientific  Research,  Boston  College,  Boston,  MA 
AFPL  F19628-95-C-01 16 

ABSTRACT: 

Seismic  moments  of  Nevada  Test  Site  (NTS)  explosions  were  determined  from  regional  surface 
wave  spectra.  Two  methods  were  used.  In  one  the  moment  is  solved  for  assuming  only  an  explosive 
source,  or  average  scalar  moment;  in  the  other  a  joint  inversion  for  an  isotropic  (explosive)  source  plus 
a  constrained  double  couple  moment  component  representing  tectonic  strain  release  (TSR).  Although 
the  general  moment  tensor  solution  to  this  joint  inversion  problem  is  non-unique,  if  some  assumptions 
are  made  concerning  the  non-isotropic  moment  components,  then  the  remaining  source  parameters  can 
be  solved  by  a  linear  least-squares  inversion  scheme.  We  examined  the  errors  in  determining  the  iso¬ 
tropic  moment  component  (A/7)  by  this  latter  method  of  constrained  linear  inversion  solutions  in  a 
canonical  study  using  a  theoretical  network  of  long-period  (6-60  sec.)  surface  wave  data.  The  network 
azimuthal  coverage  was  chosen  to  represent  that  of  a  long-period  North  American  super-network  of  55 
stations  used  for  the  actual  NTS  events.  We  compared  these  errors  in  moment  estimate  to  those 
obtained  from  surface  wave  magnitude  (M5)  and  spectral  scalar  moment  (Mq)  measurements  for  the 
same  surface  wave  observations.  For  a  ratio  of  M (eXpi)IM (eq)  less  than  1.0  we  found  that  the  inverted 
Mj  solution  is  a  much  better  estimate  of  the  actual  isotropic  moment  than  either  Ms  or  M 0,  and  the 
standard  deviation  in  this  estimate  is  substantially  less  than  that  using  the  other  two  methods  for  the 
great  majority  of  isotropic  source  +  double  couple  sources.  Even  when  the  inversion  constraints  are 
off  in  dip  and  rake  each  by  30° ,  the  mis-estimate  of  the  isotropic  moment  is  less  than  35  percent  of 
the  actual  value.  In  the  case  of  a  vertical  strike-slip  fault,  the  inverted  isotropic  moment  solution 
which  assumes  this  fault  orientation  is  exact  to  three  figures,  whereas  and  M q  under-estimate  the 
moment  by  45  percent  and  32  percent,  respectively  because  of  uneven  azimuthal  coverage. 

This  moment  tensor  inversion  method  was  applied  to  determine  the  isotropic  source  for  111  NTS 
underground  explosions  using  vertical  and  tangential  component  surface  wave  data  from  this  regional 
network.  We  also  calculated  Ms  and  M0  for  these  same  events  and  compared  the  results.  Isotropic 
source  errors  were  smallest  using  the  spectral  domain  inversion  method.  However,  this  spectral 
domain  method  cannot  attain  as  low  a  magnitude  threshold  as  the  time  domain  moment  or  A/5  method. 
The  extensive  moment  data  set  analyzed  were  combined  with  larger  yield  explosions  from  prior 
moment  studies  to  create  a  comprehensive  data  set  with  which  to  obtain  conclusive,  well-constrained 
long-period  explosion  source  scaling  relationships  at  the  separate  NTS  sub-sites. 

Regressing  on  the  results  presented  here  and  the  results  of  others  for  larger  events  with  published 
yields,  we  obtained  a  Mj  versus  yield  relation  with  which  we  were  to  estimate  the  surface  wave 
inferred  yields  of  the  1 1 1  NTS  events  before  and  after  correcting  each  event  for  it’s  sub-site  bias. 
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Preliminary  Report  on  Surface  Wave  Regionalization 

D.  G.  Harkrider  and  J.  L.  Stevens 
Institute  of  Scientific  Research,  Boston  College,  Boston,  MA 
AFPL  F19628-95-C-01 16 

SCUBED,  Division  of  Maxwell  Industries,  La  Jolla,  CA 
AFTAC  F19628-95-C-01 10 
Summary 

In  the  late  1960’s,  Kimbal  and  Kovacs  of  GEOTECH,  developed  a  lxl  degree  global  grid  of 
Rayleigh  and  Love  surface  wave  group  velocities.  The  grid  was  used  to  estimate  surface  wave  arrival 
times  for  use  in  seismic  surface  wave  association  codes;  in  other  words,  for  determining  which  of 
many  recorded  surface  waves  is  associated  with  a  particular  seismic  event.  The  emphasis  was  on  sur¬ 
face  wave  energy  in  the  17-23  second  period  band  used  for  measuring  the  surface  wave  magnitude, 
Ms ,  which  is  used  for  yield  estimation  and  as  a  discriminant  when  compared  with  the  body  wave  mag¬ 
nitude,  mb.  There  were  16  Rayleigh  wave  regions  of  which  9  were  oceans  or  shallow  seas.  South 
America,  Australia  and  Antarctica  were  considered  the  same  group  velocity  region.  Excluding  the  17- 
23  second  period  range,  the  velocities  and  more  important  the  shape  of  these  regionalized  dispersion 
curves  from  10-50  seconds  did  not  agree  with  observations  and  were  of  no  use  for  phase  matched 
filtering,  a  common  technique  for  finding  and  extracting  surface  wave  signals  from  noise. 

In  1977,  Okal  introduced  a  15x15  degree  global  grid  for  long  period  Rayleigh  waves  and  in  1989 
using  the  5x5  degree  global  tectonic  regionalization  (GTR1)  of  Jordan,  1981,  as  a  guide, 
Okal&Talandier, 1989,  presented  a  10x10  degree  grid  of  Rayleigh  group  velocities.  They  used  the 
seven  basic  regions  of  Okal  1977  (four  oceanic  regions,  shield,  trench  and  mountain)  with  velocities  in 
the  period  range  from  300  to  35  seconds.  The  ocean  velocities  were  from  Yu&Mitchell,1977,  and 
Mitchell&Yu,1980,  and  the  regionalization  was  based  on  ocean  ages.  Mitchell  and  Yu  also  determined 
separate  SV  and  SH  velocity  profiles  to  approximate  the  anisotropy  of  the  oceanic  lithosphere  and  low 
velocity  zone.  Shield,  mountain,  trench  and  subduction  zone  velocities  were  from  Nakanishi,  1981  and 
Hwang&Mitchell,  1987.  Since  1989,  Okal  and  students  have  extended  the  velocities  down  to  17 
second. 

Jordan’s  model  GTR1  has  6  regions  with  3  ocean  regions  based  on  age  and  was  presented  mainly 
as  a  geographic  framework  in  order  to  test  new  tectonic  models  or  velocity  structures  as  more  observa¬ 
tions  and  dispersion  values  become  available.  The  regional  dispersion  was  phase  velocities  from  20  to 
120  seconds  based  on  representative  curves  from  Brune&Dorman,  Landisman  et  al.  and  the  Knopoff 
UCLA  group.  Rosa  (1986)  performed  a  numerical  differentiation  of  the  scanned  phase  velocities  to 
obtain  group  velocities.  Since  our  interest  and  test  data  was  mostly  near  the  end-point  of  20 
seconds(the  least  reliable  of  the  differenced  curves),  we  were  not  able  to  use  the  majority  of  our  test 
data  source-receiver  pairs. 

Using  over  4000  source-receiver  pairs,  Stevens,  19%,  was  able  to  to  perform  a  tomagraphic 
detetrmination  of  Rayleigh  group  velocities  from  300  to  12.5  seconds  and  refine  the  Okal  dispersiom 
values  in  the  shield  and  sub-divide  the  mountain  region  into  two.  In  addition,  he  obtained  vastly 
improved  group  velocities  for  the  trench  and  subduction  region  for  a  total  of  8  distinct  dispersion 
regions.  This  10x10  degree  model  is  implemented  at  the  Prototype  International  Data  Center.  Again 
using  GTR1  as  a  guide,  we  interpolated  the  Stevens-Okal  regionalization  into  a  5x5  degree  global 
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model.  Arbitrarily  reducing  the  grid  size  for  tomographic  determined  velocities  can  create  problems 
but  in  this  case  the  inversion  path  lengths  across  the  age  boundaries  for  the  oceanic  determined  disper¬ 
sion  curves  of  Yu&Mitchell  (1979)  and  Mitchell&Yu  (1980)  were  based  on  smooth  boundaries  and  so 
in  the  oceans  the  smaller  the  grid  the  closer  to  the  actual  path  lengths  used  in  the  original  inversions. 

For  the  purpose  of  obtaining  a  more  detailed  regionalization  for  use  with  regional  dispersion  and 
magnitudes,  we  decided  to  use  Mooney’s  global  crustal  structure  model  for  calculating  shorter  period 
Rayleigh  and  Love  wave  dispersion.  Originaly  Mooney  had  a  5x5  degree  global  grid  of  52  crustal  P- 
velocity  structures  down  to  the  crust  mantle  interface(Moho).  The  crust  was  divided  into  5-7  layers.  In 
1995,  Mooney,  Laske&Masters  included  S-velocity  and  density  in  the  crustal  models.  The  result  was 
89  seismic  structures.  The  5-6  layer  P-velocities  and  interface  depths  were  obtained  from  refraction  and 
relection  profiles.  The  structures  included  Moho  velocities  and  densities.  There  can  be  a  water  or  ice 
layer  at  the  surface.  S-velorities  and  densities  were  obtained  from  surface  wave  inversions,  refraction 
profiles,  observed  relations  between  Poisson  ratio,  density,  and  depth  and  when  all  else  failed:  intuitioa 
Their  latest  model  was  1995  but  a  "final  corrected  model"  is  due  the  summer  of  1996.  Rayleigh  and 
Love  waves  were  calculated  for  periods  between  7  and  50  secs.  Because  the  structures  only  went  to 
Moho  depths,  the  oceanic  group  velocities  are  featureless  over  the  oceans  at  periods  greater  than  17 
secs.  This  resulted  from  the  common  observation  that  the  depth  of  influence  in  kms.  on  the  fundamen¬ 
tal  Rayleigh  wave  is  approximately  twice  the  period  in  secs,  and  that  all  the  ocean  structures  in  this 
model  were  identical  below  a  depth  of  20  kms.  Because  of  this  it  was  suggested  that  at  periods  near 
20  secs.,  the  Mooney  et  al.  models  be  used  for  continents  and  regions  of  thick  crust  with  a  5x5  degree 
interpolation  of  the  ocean  models  of  Okal  and  the  trenches  and  subduction  regions  of  Stevens.  This 
was  accomplished  again  using  GTR1  as  guide  and  is  our  composite  grid  of  94  dispersion  regions. 

The  6  models:  (l)Kimbal-Kovacs,  (2)Okal,  (3)Mooney,  (4)Stevens,  (5)  Stevens  5x5,  and 
(6)Mooney-Stevens  were  tested  with  the  group  arrival  times  of  the  narrow  band  filtered  period  of  max¬ 
imum  energy  in  the  range  of  17-23  secs  from  3  test  sites:  NTS,  Lop  Nor  and  Kazakh.  The  number  of 
recording  stations  ranged  from  23  to  41.  As  one  might  hope,  the  percentage  in  error  of  the  arrival 
times  was  reduced  by  each  addition  of  the  new  models.  Even  the  Mooney  model  with  its  poor  oceanic 
dispersion  was  better  than  the  Okal  model,  which  was  a  major  improvement  over  the  Kimbal-Kovacs 
model.  The  reduction  in  arrival  time  error  is  essentially  in  the  order  of  the  listed  models  except  that  the 
Mooney  model  was  better  for  the  Chinese  and  Kazakh  test-sites.  The  arrival  time  error  for  the 
Mooney-Stevens  model  was  slightly  larger  than  both  Steven’s  models  but  the  standard  deviation  of  this 
error  over  the  recording  stations  for  all  events  was  much  smaller  for  this  composite  model.  This  indi¬ 
cates  the  possibility  of  a  base  line  error  in  this  model  rather  than  in  the  regionalization  it  self.  This 
may  be  true  of  all  the  models  since  the  station  averaged  arrival  time  was  slower  than  observed  for  all 
events. 


The  Western  Wyoming  Seismic  Refraction  Profile 
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J  J  Cipar  (Phillips  Laboratory,  Earth  Sciences  Division,  Hanscom  Air 
Force  Base,  MA  01731;  617-377-3767) 
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J  A  Cipar  (R.J.  Grey  JHS,  Acton,  MA  01720) 

F  Lorenz  and  P  Lorenz  (Geoforschungszentrum,  Potsdam,  Germany) 
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E  M  O’Pray  (Analytical  Systems  Engineering  Corporation,  Burlington, 
MA  01803) 

T  Anderson  (University  of  Connecticut,  Storrs,  CT  06269) 


A  47-station  refraction  profile  was  deployed  in  western  Wyoming  to 
monitor  the  chemical  explosions  from  Project  Deep  Probe  in  August, 
1995.  The  profile  extended  from  the  edge  of  the  Wind  River  Mountains 
at  Big  Sandy,  Wyoming,  west  to  the  border  of  Wyoming  and  Idaho. 
The  first  data  set  was  collected  on  August  8th  and  9th  by  die  14  stations 
closest  to  Idaho.  The  second  data  set  was  collected  on  August  16th  and 
17th  by  the  remaining  33  stations.  Each  station  was  spaced 
approximately  3.2  km  from  each  other.  The  data  was  recorded  by  short 
period  sensors  that  were  digitized  and  stored  by  either  a  Refraction 
Technology  recorder  with  a  GPS  clock,  a  Terra  Technology  recorder 
with  WWB  radio  timing,  or  a  PDAS  recorder  with  a  GPS  clock.  All 
data  were  digitized  at  100  samples/sec  on  three  components  during  both 
periods  of  recording  with  the  exception  of  the  Terra  Technology 
recorders  which  were  digitized  at  100  samples/sec  on  the  vertical 
component. 

We  plan  to  use  the  Deep  Probe  explosion  east  of  Riverton,  Wyoming, 
recorded  by  our  profile  to  study  the  crustal  structure  beneath  the  Wind 
River  Mountains  in  south-western  Wyoming  as  well  as  the  Wyoming 
and  Salt  River  Ranges  near  Idaho.  The  Deep  Probe  explosions  in 
Canada  and  the  US  will  be  used  to  explore  the  deep  structure  of  western 
North  America. 
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The  Wyoming  Large- Aperture  Seismic  Array 

John  Cipar  (Earth  Sciences  Division,  Phillips  Laboratory,  Hanscom 
AFB,  MA  01731;  617-377-3767;  e-mail:  cipar@doc.plh.af.mil) 
Stephan  Koester  (Boston  College,  Institute  for  Scientific  Research; 
617-377-2652;  e-mail:  koester@doc.plh.af.mil) 

James  Steeves,  Eileen  O'Pray  (Analytical  Systems 
Engineering  Corporation,  Burlington,  MA  01803) 

Thomas  Anderson  (University  of  Connecticut,  Storrs,  CT); 

James  Cipar  (R  J.  Grey  JHS,  Acton,  MA) 

A  five  station  seismic  array  has  been  installed  in  southwestern 
Wyoming  since  late  July,  1995.  The  array  is  located-in  the  Green 
River  basin,  west  of  the  Wind  River  Mountains.  The  five  stations 
are  arranged  in  the  form  of  a  cross  with  the  central  station  near  Big 
Piney,  Wyoming.  Each  outlying  station  is  located  about  50  km  from 
the  central  site.  Broadband  and/or  long-period  sensors  are  installed 
at  each  station  on  a  concrete  pad  surrounded  by  an  insulated, 
wooden  enclosure.  Data  are  recorded  on  Refraction  Technology  data 
acquisition  systems  equipped  with  530-Mbyte  disks  and  GPS  clocks. 
For  the  first  two  weeks  of  recording,  during  the  period  of  the  - 
Canadian/US  Deep  Probe  experiment,  data  were  recorded 
continuously  at  40  samples/sec.  Presently  the  array  is  recording 
continuously  at  20  sps  (10  sps  at  2  sites).  Power  is  supplied  either 
from  solar  panels  or  by  line  power  at  certain  sites.  Current  plans  are 
to  operate  the  array  into  early  October,  1995. 

We  plan  to  use  the  data  from  this  array  to  study: 

-  regional  wave  propagation  in  the  central  Rocky  Mountains 

-  spatial  variability  of  earthquake/explosion  discriminants 

-  crustal  structure  beneath  the  array 

Examples  of  data  from  the  array  such  as  the  August,  1995  Chinese 
nuclear  test,  regional  earthquakes,  and  Deep  Probe  and  mining 
explosions  will  be  displayed. 
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ABSTRACT 

In  June  of  1995,  the  Earth  Sciences  Division  of  the  Air  Force  Phillips  Lab,  with  survey  equipment 
from  the  University  of  Delaware  and  assisted  by  the  Kansas  Geological  Survey  and  Elohi  Geophysics, 
conducted  a  geophysical  site  characterization  of  the  SERDP-funded  Groundwater  Remediation  Field 
Lab  (GRFL)  located  at  Dover  AFB,  Delaware  and  administered  by  Applied  Research  Associates  for 
USAF  Armstrong  Lab.  Seismic  data  were  collected  in  order  to  1)  compare  the  results  using  three 
different  compressional  sources  and  2)  cover  the  field  site  well  enough  to  characterize  the  seismic 
response  of  the  shallow  subsurface.  This  paper  will  focus  primarily  on  the  first  of  these  two  goals. 

Seismic  data  were  collected  along  three  north-south  profiles  set  10  meters  apart,  each  profile  with  a 
different  compressional  source:  a  5.5kg  sledgehammer,  a  12-gauge  firing  rod  from  Betsy  Seisgun  Inc. 
shooting  150  grain  blanks,  and  a  portable  piezoelectrically  driven  vibrator,  developed  by  Elohi 
Geophysics,  operating  with  a  90Hz-450Hz  sweep.  An  east-west  cross  line  was  collected  using  the 
sledgehammer  source  in  order  to  tie  the  three  profiles  together.  A  laser  theodolite  provided  station 
location  and  elevation  control.  The  primary  targets  were  the  water  table  (that  had  been  marked  on 
maps  at  a  depth  of  about  3  meters)  and  a  sand-clay  interface  at  about  15  meters  depth.  We  collected 
24-channel  CMP  data  using  a  half  meter  spacing  of  both  source  and  100Hz  geophones.  An  end-on 
spread  geometry  was  used,  with  a  1  meter  offset  between  source  and  nearest  geophone.  Field  QC 
after  initial  walkaway  noise  testing  with  each  source  did  not  show  any  one  source  to  be  outstanding 

We  have  associated  the  strongest  reflection  event  with  the  sand-clay  interface  interpreted  as  the  top 
of  the  aquitard.  A  practical  early  result  of  the  seismic  survey  showed  the  water  table  to  be  at  over  8 
-  meters.  Seismic  data  comparison  in  this  study  is  based  on  spectral  content,  total  energy  and  signal- 
to-noise  ratios,  as  well  as  a  discussion  of  coherency  of  the  primary  reflection  event  at  the  water 
table.  The  problem  with  the  water  table  being  deeper  than  expected  is  that  the  water  table  reflection 
may  interfere  with  the  primary  seismic  target,  the  sand-clay  interface.  With  a  wavelength  of  about  4 
meters  at  100  Hz,  interpretation  of  the  data  must  take  into  account  the  possible  interference  of  the 
two  reflections  in  the  seismic  images. 

With  a  surface  velocity  of  400m/sec,  the  first  Fresnel  zone  for  100Hz  signals  at  15  meters  depth  is 
about  5.5  meters  under  each  seismic  line,  therefore  overlapping  between  profiles.  Thus,  despite  the 
separation  of  the  three  lines,  they  are  sampling  similar  regions  of  the  target  area.  Nevertheless, 
initial  inspection  of  the  seismic  shot  gather  data  showed  that  they  are  characterized  by  rapid 
variations  in  amplitude  and  phase  across  short  distances. 

Both  the  firing  rod  and  the  Vibratory  source  gave  an  initial  look  at  the  near  surface  during  data 
acquisition  via  the  use  of  augers  necessary  for  deployment  of  these  sources.  The  site  had  rapid  lateral 
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changes  in  the  upper  meter.  Clays  and  gravel  stringers  with  lateral  variability  on  the  order  of  a  half 
meter  were  the  norm.  Cone  penetrometer  data  suggest  that  this  heterogeneity  extends  deeper  as 
well.  Our  comparison  of  the  data  acquired  using  the  different  sources  is  with  the  caveat  that  there  is 
an  extreme  variability  in  the  near  surface. 


BACKGROUND  AND  SITE  CHARACTERISTICS 

Dover  Air  Force  Base,  Dover,  Delaware,  is  the  site  of  the  Groundwater  Remediation  Field  Laboratory 
(GRFL),  a  national  test  site  administered  by  the  Air  Force  Armstrong  Laboratory  (AL)  with  funds 
from  the  Strategic  Environmental  Research  and  Development  Program  (SERDP),  a  partnership 
between  DOD,  DOE  and  EPA.  In  1995,  AL  and  their  contractor  Applied  Research  Associates  (ARA) 
had  laid  the  groundwork  for  the  GRFL  which  will  provide  an  infrastructure  in  which  a  number  of  test 
cells  will  be  able  to  operate  simultaneously.  The  GRFL  is  intended  to  offer  a  place  where  researchers 
may  test  remediation  technologies  that  address  a  variety  of  LNAPL  and  DNAPL  targets  in  the 
subsurface.  Geophysical  characterization  of  the  site  was  begun  by  The  Air  Force  Phillips  Lab  (PL) 
and  the  University  of  Delaware  in  March- June  1995,  collecting  ground  penetrating  radar  (GPR)  at 
multiple  source  frequencies,  shallow  high  resolution  seismic  reflection  data  and  terrain  conductivity. 
ARA  continued  through  the  end  of  the  summer  collecting  soil  samples,  cone  penetrometer  data,  and 
more  GPR. 

The  test  site  consists  of  an  unconfined  aquifer  extending  to  a  depth  of  15-18  meters,  with  a  clay 
aquitard  below.  The  water  table  is  at  about  8.5  meters.  The  surface  of  the  site  is  a  3.5  acre  grassy 
field  that  has  a  gentle  dip  to  the  east,  with  only  minor  topographic  undulations  of  less  than  .5  meter. 
Initial  goals  of  the  seismic  survey  were  to  characterize  the  seismic  response  of  the  site  and  to  collect 
reflection  profile  data  coincident  with  specific  lines  of  the  more  densely  spaced  GPR  data.  Previous 
work  in  coastal  plain  sedimentary  environments  (Miller,  et  al.,  1986)  suggested  that  most 
compressional  source  types  should  be  adequate  for  the  seismic  reflection  data  acquisition,  and  our 
initial  plan  was  to  choose  the  best  source  for  making  a  comparison  between  the  seismic  and  the  radar 
data.  After  preliminary  source  parameter  testing  in  April  1995,  we  determined  that  a  seismic  source 
comparison  at  this  site  might  be  worthwhile.  The  seismic  reflection  data  were  acquired  in  June  1995. 


DATA  ACQUISITION  AND  FIELD  QC 

All  seismic  data  were  collected  with  a  24-channel  Geometries  StrataView  recorder.  The  geophones 
were  single  component,  100Hz  Mark  Products  L40A  instruments  with  14cm  spikes.  Standard  CDP 
cable  and  roll  switch  were  used  to  collect  the  reflection  data.  Figure  1  shows  the  field  site  geometry 
and  layout  of  CMP  lines.  We  collected  24-channel  CMP  data  using  a  half  meter  spacing  of  both 
Source  and  geophones.  An  end-on  spread  geometry  was  used,  with  a  1  meter  offset  ^between  source 
and  nearest  receiver.  Three  compressional  sources  were  tested  and  used  for  CMP  profiling:  a  5.5kg 
sledgehammer,  a  Betsy  Seisgun,  Inc.  12-gauge  firing  rod,  and  the  Earth  Reaction  Seismic  Source 
(ERSS)  portable  vibrator  from  Elohi  Geophysics. 

In  one  day  we  collected  332  source  points  of  sledgehammer  data  along  a  line  through  the  central 
portion  of  the  field.  At  each  source  point  we  stacked  eight  hits  of  the  hammer,  for  a  total  of  over 
2600  blows.  Figure  2  shows  a  representative  shot  gather  from  the  sledgehammer  CMP  collection. 
The  air  wave  and  a  slow  direct  wave  (~250m/sec)  overlap  with  the  ground  roll  at  near  offsets.  A 
shallow  refractor  at  ~475m/sec  is  strong  in  the  far  offset  traces.  Although  source  repeatability  was 
not  optimum  due  to  four  different  hammer  operators  and  the  length  of  day  in  the  field,  stacking 
eight  blows  per  shot  point  did  a  very  good  job  of  equalizing  source  wavelet  and  signal  strength. 

Ten  meters  to  the  east  of  the  sledgehammer  line  we  collected  187  shot  points  with  a  12-gauge  firing 
rod  shooting  150  grain  blank  loads.  Holes  were  drilled  about  one  meter  deep,  and  three  different 
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packing  techniques  were  tried:  no  packing  (air-filled),  water-filled  and  dirt-filled.  There  was  not 
much  difference  noticeable  in  the  shot  data  between  water-filled  holes  and  dirt-filled  holes.  The  bulk 
of  the  CMP  data  were  collected  with  dirt-filled  holes  for  source  coupling.  Figure  3  shows  a 
representative  shot  gather.  Whereas  the  sledgehammer  had  a  large  ratio  of  ground  roll  to  body  wave 
energy,  the  firing  rod’s  strongest  arrival  is  the  475m/sec  refraction  event  that  on  some  records  is 
even  stronger  than  the  air  wave. 

Ten  meters  to  the  west  of  the  sledgehammer  line  we  collected  228  source  points  of  CMP  data  with 
the  ERSS  vibrator  (Figure  4).  The  ERSS  was  operated  with  a  linear  sweep  90Hz  to  450Hz  three 
seconds  long.  Each  source  point  was  a  stack  of  eight  sweeps.  Deployment  of  the  ERSS  is  by 
hydraulic  augering  in  ground  anchors  and  hooking  the  source  to  them  one  at  a  time.  The  hydraulic 
auger  also  is  used  to  remove  anchors.  The  ERSS  weighs  about  90  pounds  and  the  hook  clamping 
system  is  set  to  generate  300  pounds  of  hold-down  force.  At  test  time,  a  pilot  was  not  recorded  with 
each  sweep  so  uncorrelated  data  were  recorded  and  cross  correlation  done  with  a  synthetic  sweep 
during  processing.  The  FK-spectra  of  the  vibrator  records  show  low  groundroll  and  airwave  energy 
relative  to  that  of  body  waves.  The  shot  records  (Figure  5),  as  with  most  all  vibrator  data,  do  not 
show  clear  refraction  arrivals. 

Planting  the  geophones  was  not  difficult  over  all;  however  deployment  of  the  ERSS  and  the  Betsy 
Seisgun  required  penetrating  deeper  into  subsurface.  As  evidenced  by  the  auger  work  for  each  of  these 
sources,  the  variability  of  the  upper  meter  of  the  subsurface  was  extreme  on  a  scale  of  a  half  meter 
(from  one  shot  point  to  the  next  shot  point).  Clay  stringers  and  gravel  lenses  were  the  norm  across 
the  entire  field,  and  cone  penetrometer  data  indicate  this  variability  extends  well  below  the  upper  half 
meter.  The  lateral  variability  makes  it  difficult  to  compare  directly  the  stacked  sections  from  each 
source  relative  to  strength  and  continuity  of  common  reflection  events.  Nevertheless,  we  will  discuss 
these  issues  as  well  as  make  some  observations  regarding  the  resolution  of  reflectors  in  the  data. 


DATA  PROCESSING  AND  INTERPRETATION 

Each  CMP  data  set  for  this  comparison  went  through  the  following  processing  history: 

Bandpass  filter  (zero-phase  Butterworth  80-300Hz  with  Hamming  taper) 

F-K  filter  to  remove  airwave  and  groundroll 

First  arrival  mute  (direct  wave  and  shallow  refraction) 

NMO  correction  (after  constant  velocity  stacks  to  get  best  velocity  functions) 

CMP  stack 

Signal  to  noise  enhancement  with  single-adjacent  trace  mixing 
Display  with  50ms  AGC 

Figure  6  displays  the  three  CMP  profiles  for  comparison.  The  sledgehammer  and' vibrator  lines  are 
cut  short  to  coincide  with  the  length  of  the  firing  rod  line.  In  the  sledgehammer  CMP  profile  (Figure 
6a),  there  are  three  major  reflection  events  at  ~34-38ms,  ~40-48ms  and  at  ~25-28ms.  In  the  south 
end  of  the  line,  the  deeper  reflections  dominate,  while  at  about  CMP  location  350  the  shallower 
event  becomes  prominent.  Both  of  the  deeper  events  shallow  somewhat  from  south  to  north  as  they 
approach  the  middle  of  the  profile,  and  the  shallower  event  has  a  similar  dip. 

The  shallow  refractor  witnessed  in  the  shot  gather  data  at  about  lm  depth  corresponds  with  what 
might  be  the  bottom  of  a  tilled  zone  for  this  field  that  at  one  time  had  been  farmland.  The  ~36ms 
reflection  event  corresponds  to  a  reflection  off  the  water  table  at  8.1m.  The  ~44ms  reflection  event 
corresponds  to  the  top  of  the  clay  confining  layer  at  14m  depth.  Data  processing  of  all  three  seismic 
lines  suggested  that  the  near  surface  velocity  increased  toward  the  north.  During  data  acquisition  it 
was  noted  that  the  north  end  of  the  field  was  extremely  different  at  the  shallowest  level.  Geophones 
and  survey  markers  were  harder  to  plant  due  to  much  more  compacted  surface.  The  different  nature 
of  the  shallowest  layering  might  be  due  to  the  construction  of  the  running  track  to  the  northeast  of 
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the  study  area  (Figure  1).  This  velocity  change  from  south  to  north  could  explain  the  apparent 
shallowing  of  the  reflection  events  that  are  more  likely  horizontal  or  even  deepening  (as  limited 
cone  penetrometer  data  suggests). 

The  firing  rod  CMP  profile  (Figure  6b)  mirrors  the  sledgehammer  line,  with  the  42-44ms  reflection 
event  dominating  the  record.  The  profile  stops  just  about  at  the  point  where  the  shallower  reflection 
event  becomes  prominent  in  the  sledgehammer  line.  North  of  CMP  location  245  the  firing  rod  line 
shows  this  shallower  event  as  well.  The  top-of-clay  reflection  is  much  more  continuous  across  the 
profile,  with  better  signal  to  noise  evident  than  in  the  sledgehammer  line.  The  water  table  reflection 
is  just  as  discontinuous  as  in  the  sledgehammer  profile.  There  is  a  suggestion  in  the  firing  rod  profile 
of  something  deeper  than  the  top-of-clay;  however,  it  seems  to  track  the  prominent  event  above  it, 
and  thus  is  very  likely  a  side  lobe  of  the  top-of-clay  event. 

With  low  groundroll  and  airwave,  and  little  direct  wave,  the  vibrator  data  offers  the  possibility  for  the 
shallowest  imaging,  although  the  first  reflector  is  not  seen  until  two-way  time  of  ~25ms  in  the  north 
portion  of  the  profile  (Figure  6c)  as  seen  in  the  sledgehammer  and  firing  rod  profiles.  The  vibrator 
profile  shows  the  water  table  reflection  and  the  top-of-clay  as  with  the  other  sources,  however  the 
water  table  reflection  appears  to  be  the  stronger  and  more  continuous  of  the  two. 

Figure  7  shows  amplitude  spectra  for  each  stacked  section  using  a  60  ms  window,  starting  at  5ms  two- 
way  time  for  the  signal.  Comparing  the  average  spectra  after  stack  across  each  section,  the  vibrator 
data  is  the  most  broadband  with  a  peak  at  125Hz.  The  sledgehammer  offers  the  next  best  frequency 
content,  with  a  peak  at  93Hz  and  with  more  usable  information  out  to  200Hz  than  with  the  firing 
rod.  The  firing  rod  peaks  at  85Hz,  and  rolls  off  a  little  more  quickly  at  higher  frequency  than  does 
that  of  the  sledgehammer. 

Presumably,  the  lower  resolution  in  the  impulsive  data  records  result  in  a  little  more  interference 
between  die  two  primary  reflection  events.  This  may  explain  the  difference  between  the  impulsive 
sources  that  show  the  top-of-clay  event  as  the  most  coherent  and  the  vibrator  profile  that  shows  the 
water  table  as  the  most  coherent  event. 

Amplitude  spectra  for  single  traces  in  each  stacked  section  vary  quite  a  bit.  The  three  sources  have 
similar  maximum  amplitudes  on  a  trace-by  trace  basis  (Figure  7).  The  sledgehammer  was  the 
strongest  after  summing  eight  blows,  whereas  the  vibrator  was  the  weakest.  Average  amplitude 
spectra  across  each  section  show  the  firing  rod  has  a  more  consistent  source  spectrum,  retaining  its 
amplitude  and  leaving  an  average  amplitude  that  is  about  twice  the  strength  of  the  sledgehammer  and 
four  times  the  strength  of  the  average  for  the  vibrator.  Using  a  time  window  deeper  in  the  profiles 
(90ms- 150ms)  to  get  an  estimate  of  noise  characteristics,  the  average  noise  amplitude  peaks  at 
1157Hz  for  the  sledgehammer  and  the  firing  rod,  and  at  202Hz  for  the  vibrator.  Using  the  peak 
values  from  the  averaged  spectra  in  Figure  7  to  get  the  best  signal  and  the  worst  noise,  the  signal  to 
noise  ratio  (S/N)  for  the  firing  rod  is  the  greatest  at  4.3.  The  sledgehammer  S/N=2.1  and  the  vibrator 
S/N=1.7  are  close  to  each  other  and  at  less  than  half  that  of  the  firing  rod. 


CONCLUSIONS 

Although  this  work  constitutes  a  case  study  for  an  area  where  it  is  somewhat  difficult  to  achieve 
superb  seismic  results,  it  also  offers  a  good  comparison  of  three  near-surface  seismic  sources  for  CMP 
data  collection  and  interpretation.  In  this  paper  we  have  focused  specifically  on  the  source 
comparison.  Future  reports  will  cover  the  more  detailed  interpretation  of  the  data  and  the  seismic 
characterization  of  the  Dover  AFB  site. 

With  respect  to  resolution  of  reflection  events  in  this  study,  the  most  broadband  signal  was  obtained 
by  the  portable  vibratory  source.  The  vibrator  offered  continuity  of  reflection  events  similar  to  that 
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of  the  impulsive  seismic  sources,  although  events  were  difficult  to  interpret  from  shot  gathers  and  no 
acquisition  QC  could  be  done  directly  on  the  recorder  while  collecting  the  uncorrelated  data.  The 
firing  rod  source  offered  by  far  the  best  signal  to  noise  ratio  and  arguably  the  best  lateral  continuity 
of  reflection  events.  The  eight  blows  of  the  sledgehammer  offered  a  source  strength  that  was  in 
between  those  of  the  firing  rod  and  the  vibrator,  with  a  frequency  content  somewhat  better  than  the 
firing  rod.  The  sledgehammer  was  certainly  the  fastest  source  for  collecting  the  most  CMP  data  in  a 
limited  time,  and  was  no  more  labor  intensive  than  the  other  two  sources  that  required  the  use  of 
two-person  augers.  As  each  source  excelled  in  one  certain  aspect,  there  was  no  single  winner  in  our 
comparison.  Specific  survey  and  site  limitations,  and  desired  targets  and  results,  must  still  dominate 
the  choice  of  compressional  source  for  near-surface  exploration. 
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DOVER  AIR  FORCE  BASE  TEST  SITE  _ JUNE  1995 


Figure  1 .  Geophysical  survey  layout  at  Dover  AFB 
Groundwater  Remediation  Field  Lab  site 
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Figure  4.  Deploying  the  Earth  Reaction  Seismic  Source  (ERSS) 
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Figure  5.  Example  shot  record  for  the  portable  vibrator.  Filter  parameters  as  in  Figure  2 


Time  (sec) 


CMP  position 


North 


25  50  75  100  125  150  175  200  225  250  275  300  325  350 


0.0 


0.1 


25  50  75  100  125  150  175  200  225  250  275  300  325  350 


Figure  6.  CMP  profiles  for  (a)  sledgehammer,  (b)  firing  rod  and  (c)  vibrator 
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1.  INTRODUCTION 


The  idea  of  cellular  automata  dates  back  to  the  work  of  von  Neumann  and  Ulam  in  the 
1940s  (von  Neumann  1966)  who  were  looking  for  simple  rules  of  spatial  and  temporal  evolution 
to  mimic  collective,  complex  behavior  of  biological  systems. 

The  first  goal  was  to  provide  a  theory  for  how  an  artificial  life,  capable  of  self¬ 
reproduction  could  be  constructed,  starting  with  the  simplest  rules  possible.  So  appeared  the  idea 
of  a  collection  of  very  simple  finite-state  machines  (cellular  automata)  with,  for  simplicity,  only 
binary  values,  the  state  of  each  depending  only  of  the  state  of  its  immediate  environment 

The  natural  space  on  which  to  put  all  this  is  a  lattice  with  elementary  finite-state  machines 
placed  at  the  vertices.  The  rules  for  updating  this  array  of  small  machines  (lattice  cellular 
automata)  can  be  done  concurrently  in  one  clock  step,  that  is,  in  parallel. 

In  the  last  decades  the  interest  from  reproducing  biological-like  behavior  automaton  has 
shifted  towards  physics  and  computation:  simulation  of  partial  differential  equations  (Zuse  1970), 
time-reversible  automata  (Margolus  1984),  simulation  of  quantum-mechanical  phenomena 
(Feynman  1982),  “statistical  mechanics”  of  cellular  automata  (Wolfram  1983),  cellular-automata 
as  discrete  dynamical  systems  (Vichniac  1984)  and  alternative  to  partial  differential  equations 
(Toffoli  1984). 

Special  cases  of  cellular  automata  are  capable  to  describe  complex  collective  behavior  and 
the  simulation  of  fluid  dynamics  is  one  in  which  the  cellular  automata  proved  they  are 
complementary  and  powerful  tools  to  model  phenomena  that  would  normally  be  the  exclusive 
domain  of  partial  differential  equations. 

This  class  was  named  lattice-gas  cellular  automata  and  have  not  only  become  a  “toy 
models”  for  the  exploration  of  the  microscopic  basis  of  hydrodynamics,  but  also  tools  for  the 
numerical  study  of  complex  problems  in  fluid  mechanics. 

While  the  lattice  gas  may,  in  principle,  be  used  for  nearly  any  problem  in  hydrodynamic 
simulation,  many  of  the  most  successful  applications  have  involved  either  a  complex  fluid,  a 
complex  geometry,  or  both. 

In  the  following  we  first  introduce  the  FHP  model  (Boon  1991,  Clouquer  et  al.  1987, 
Comubert  et  al.  1991,  Dubrulle  et  al.  1991,  Frisch  and  Rivet  1986,  Frisch  et  al.  1986,  Hayot 
1987,  Henon  1987a,  1987b,  1992,  d’Humieres  and  Lallemand  1987,  1990,  Kadanoff  et  al.  1987, 
Orszag  and  Yakhot  1986,  Rivet  and  Frisch  1987,  Rivet  et  al.  1987,  Shimomura  et  al.  1987, 
Sommers  and  Rem  1991,  Vichniac  1990,  Zanetti  1990,  1991),  the  first  of  the  wide  class  of  fluid 
models  known  as  lattice-gas  automata  (LGA),  describe  its  hydrodynamic  limit,  and  illustrate  its 
ability  to  simulate  the  Navier-Stokes  equations  for  both  2D  and  3D  flow.  Further  on,  we  will 
introduce  the  later  developments  of  the  LGA’s,  namely,  the  lattice  Boltzmann/LGB  (Benzi  et  al. 
1992,  Frisch  1991,  Higuera  et  al.  1989,  Higuera  and  Jimenez  1989,  Higuera  and  Succi  1989, 


29 


d’Humieres  1992,  McNamara  and  Zanetti  1988,  Rivet  and  Frisch  1988),  and  Bhatnagar-Gross- 
Krook  models  /LBGK  (Alexander  et  al.  1992,  Behrend  et  al.  1994,  Boghosian  and  Taylor  1995, 
Chen  et  aL  1992,  Chen  et  al.  1994,  Qian  et  al.  1992,  Qian  and  Orszag  1993). 

In  the  remainder  of  this  overview,  we  present  some  of  the  main  applications  and 
numerical  experiments  that  has  been  performed  with  lattice-gas  models:  binary  fluid  mixture  and 
phase-separation  transitions  (Alexander  et  al.  1992,  Appert  1993,  Appert  et  al.  1991,  Burges  and 
Zaleski  1987,  Chan  and  Liang  1990,  Rothman  and  Keller  1988,  Rothman  1992,  Rothman  and 
Kadanoff  1994,  Rothman  and  Zaleski  1994),  multi-phase  flow  through  porous  media  (Behrend 
1995,  Ferreol  and  Rothman  1995,  Flekkoy  and  al.  1992,  Gunstensen  1992,  Gunstensen  and 
Rothman  1992,  Gutfraid  and  Ippolito  1995,  Lutzko  et  al.  1990,  Olson  1995,  Pot  1994,  Rothman 
1988  and  1990,  Shan  and  Chen  1992,  Schwartzer  1995),  diffusion  (Burges  and  Zaleski  1987, 
Brieger  and  Bonomi  1991,  Chopard  and  Droz  1988,  d’Humieres  et  al.  1988,  van  der  Hoef  and 
Frenkel  1991,  McNamara  1990,  Qian  et  al  1992),  reaction-diffusion  equation  (Dab  et  al.  1990, 
Liu  and  Goldenfeld  1991,  Wells  et  al.  1991,  Gerits  and  Ernst  1993,  Lawniczack  and  al.  1991), 
magneto-hydrodynamics  (Chen  et  al.  1992,  Hatori  and  Montgomery  1987,  Martinez  et  al.  1993), 
free  boundary  flows  (Cliffe  et  al.  1991,  Mujica-Femandez  1991,  Eggels  1995,  Eggels  and 
Sommers  1995,  Skordos  1995,  Teixeira  1992,  Alexander  et  al.  1992),  transport  of  particle 
suspensions  (Behrend  et  al.  1994,  Behrend  1995,  Ladd  1990,  Ladd  1991,  Ladd  1994a,  Ladd 
1994b,  Ladd  et  al.  1995). 

2.  Lattice-gas  cellular  automaton  models  of  fluids 

Fluid  dynamics  is  an  especially  good  domain  for  a  cellular  automaton  formulation 
because,  in  this  case,  at  the  microscopic  level,  we  have  many  simple  atomic  elements  colliding 
rapidly  with  simple  interactions  which  coincides  with  our  intuitive  picture  of  dynamics  on  a 
cellular  space.  While  at  the  microscopic  level,  physical  fluids  consists  of  discrete  particles,  on  a 
large  scale,  they  seem  continuous  and  can  be  described  by  the  partial  differential  equations  of 
hydrodynamics,  and  in  fact,  the  form  of  these  equations  are  quite  insensitive  to  microscopic 
details  (changes  in  molecular  interactions  laws  can  affect  the  transport  parameters,  such  as 
viscosity,  but  do  not  alter  the  basic  form  of  the  macroscopic  equations). 

This  remarkable  similitude  has  not  only  had  implications  for  statistical  mechanics  and 
kinetic  theory  but  also  for  the  numerical  simulation  of  certain  hydrodynamic  flows. 

2. 1.  Fluid  dynamics 

It  will  be  very  useful  later,  in  the  construction  of  the  discrete  model  of  fluids,  to  analyze 
the  complementarity  between  the  microscopic  (kinetic)  approach  and  the  continuum 
(macroscopic)  approach  in  the  recovery  of  the  main  equations  of  fluid  dynamics. 
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2. 1. 1.  The  continuum  approach 


The  usual  way  of  deriving  the  Navier-Stokes  equations  (Batchelor  1967,  Landau  and 
Lifschitz  1959,  Reif  1988)  is  to  start  from  the  conservation’s  laws  of  mass  and  momentum, 
applied  to  small  (but  big  enough  to  avoid  their  microscopic  structure)  surfaces  and  volumes  in  a 
fluid  field. 

If  v(x,  t)  is  the  vector  field  of  the  “macroscopic”  fluid  cells  and  X  a  generalized  volume  in 
this  field,  then  the  total  mass,  flowing  out  from  X  has  to  balance  the  fluid  flow  through  the 
surface  S  (Wolfram  1986,  Hasslacher  1987,  Rothman  and  Zaleski  1996)  r 

i  pvdS  =  [  pvndS  (1) 

JdL  JBL 

where  v  is  the  velocity  vector,  5Zis  the  boundary  of  X,  S  is  the  surface  and  n  is  the  normal 
vector  of  the  surface  (positive  in  the  outward  direction). 

From  the  Gauss  law  (<j>^A-dS  =  J^V  A)dV  where  V  is  a  volume  in  R3  and  A  a  n- 

differential  form)  we  have: 

pv  dS  =  —  a  f  pv  n dS  =  [  v  •  (pv)  dV  (2) 

J  <9X  lJdL  J2 


From  which  we  obtain  the  continuity  equation  (where  9t  stands  for  — ): 

at 


dtp  +  V(pv)  =  0 


(3) 


Now,  if  p  is  the  pressure  exerted  by  the  fluid  on  the  unit  surface  area  of  an  enclosed 
dv 

volume,  and  F  =  m  a  =  pV — ,  is  the  Newton’s  second  law,  the  total  force  acting  on  a  volume  of 

dt 

fluid  due  to  the  remainder  of  the  fluid  can  be  written  as:  F  =  -  ®  p  dS. 

Using  again  the  Gauss  law,  we  obtain: 

F  P  ^  =  (4) 


and  therefore  the  force  acting  on  the  unit  mass  of  fluid  is  (-Vp). 

Then,  from  Newton’s  second  law  for  the  unit  mass  of  fluid,  and  from  (4)  we  have: 


*•=-?/>=  p£  =  p{<?,v+(v-V)v} 


(5) 


31 


where  we  used  that  dv/dt  is  a  total  derivative  (e.g.  dv  =  ^-dt  +  ^dx+^-dy+  — dz). 

at  ax  ay  dz 

Finally,  from  Eqs.  (4)  and  (5)  we  obtain  the  Euler's  equation  for  an  ideal  dissipation-free 
fluid,  as: 


d,v  +  (v-  V)v  =  Vp  (6) 

P 

Equations  (3)  and  (6)  are  a  set  of  inviscid  equations  for  an  ideal  compressible  flow,  the 
pressure  variation  law  being  given  by  the  thermodynamic  equation  of  state. 

For  real  fluids  we  have  to  take  into  account  the  dissipative  effect  of  the  viscosity  and 
therefore,  for  this  purpose,  it  is  useful  to  analyze  the  momentum  flux  time  change  rate. 

If  pv  is  the  momentum  of  the  fluid  passing  through  the  elementary  volume  dV,  then  its 
time  rate  of  change  can  be  expressed  as: 

dt(pVi)  =  {dtp)-Vi+p(dtVi)  (7) 

where  vi  are  the  velocity  components  on  the  space  dimensions. 

From  (3),  (6)  and  (7)  we  have  then: 

dt  (P  vi )  =  “  di  P  ~  P  vjfc  •  dk  vi  ~  vi  •  dk  (pv *)  (8) 

where  we  used  the  Einstein  summation  convention  (repeated  subscripts  imply  a  summation). 
Equation  (8)  is  equivalent  with: 

dt{pvi)  =  -diP-dk(PviVk )  (9) 

which  can  be  written  as: 

dt{pvi)  =  ~dk{pSik+pvivk)  (10) 

where  8%  is  the  Kronecker  delta,  and  the  expression  II&  =  P  8ik  +  P  vi  vk  is  the  inviscid 
momentum  flux  density  tensor  (Landau  and  Lifshitz  1959). 

Finally,  we  can  write  the  momentum  balance  (Euler's  equation)  in  the  condensed  form: 


IK) 


(ll) 


where  II/*  is  35  above  and  has  the  significance  of  the  i*  component  of  momentum  flowing  in  the 
direction. 

We  can  generalize  this  expression  for  the  case  of  viscous  flows  as: 

*&>».•) —Atftla)  (12) 

where  11/*  is  the  momentum  stress  tensor,  which  accounts  also  for  the  viscous  effects  in  the 
fluid. 

Usually,  II/*  is  written  in  the  form: 

ik-ni+ng*  (B) 

where  n&  has  the  same  meaning  as  before  and  1"I]*5C  is  the  viscous  stress  tensor  which  accounts 
for  the  viscosity  effects  (shear  and  compression)  in  a  Newtonian  fluid.  To  describe  the  effects  of 
viscous  stress,  one  may  introduce  an  unknown  tensor  ft/*,  so  as  II/*5C  =  -p^/*,  and  therefore, 
the  momentum  stress  tensor  now  reads: 


II;*  =  II/*  -  P  1tik  =  P  5ik  +  p  V/  V*  -  P  ft/*  =  ft/*  +  p  V/  V*  (14) 

where  ft/*  is  the  viscosity  stress  tensor  and  ft/*  =  pSik  ~P  Xik  is  called  the  stress  tensor. 

The  general  form  of  ft/*  can  be  deduced  from  some  basic  facts:  in  the  assumption  of 
small  velocity  gradients  «  <?*v/,  and  ft/*  =  0  for  v  =  0  and  under  rotation  (uniform  rotation 

doesn’t  produce  transport  of  momentum  overall). 

The  general  form  that  fulfill  these  assumptions  can  be  written  as: 


n\k  =  A{dkVi  +  divk)  +  B 8ik  dj  v j 


(15) 


and,  therefore,  the  viscous  stress  tensor  will  be: 

UWC  =  ~PA(dk  v/  +  di  v*)  -  pBSik  djVj  =  -p(<9*v/  +  <9/v*)-£  SucdjVj  (16) 

where  /i  is  the  dynamic  shear  viscosity  and  £  is  the  bulk  viscosity. 

Now,  from  Eqs.  (12),  (13),  and  (16)  we  obtain: 

dt  {Pvi)  =  ~d{[p  5ik  +  P  Vi  v*  -  p{dic  v/  +  <9/v*)  -%5ikd  jVj\  (17) 
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that  finally  can  be  written  as: 


dt  (Pvi) + dk  (p  v;  vjt)  -~diP+dk  [p{dk  v» +  <9;  vjfc)] +  di(%djVj)  (1 8) 

which  is  the  momentum  conservation  equation  for  a  compressible  viscous  fluid. 

In  the  assumption  of  an  incompressible  fluid  ( p=p0 ,  constant)  and  <9*  v*  =  0 ,  the  Navier- 

Stokes  equations  will  read: 

dt  Vi+ (v  j  d  j)  v,-  =  “  di  P  +  v^i  (19) 

where  v  =  —  is  the  kinematic  shear  viscosity.  In  vector  form,  Eq.  (19)  will  be: 

P 

dtv+(v-V)v  =  ~'Vp  +  vy2v  (20) 

2. 1. 2.  The  kinetic  approach 

The  Navier-Stokes  equations  of  fluid  dynamics,  are  highly  non-linear  and  therefore  they 
allows  for  analytic  solutions  only  for  special  cases.  Usually  they  can  be  resolved  only  by 
approximation  techniques,  the  analysis  of  Navier-Stokes  equations  still  remaining  a  mater  of 
ingenious  techniques  covering  only  particular  regimes  and  specific  geometry  (Frisch  et  al.  1987, 
Hasslacher  1987,  Wolfram  1986).  Therefore,  complementary  technique  had  to  be  developed,  in 
order  to  approach  other  real-world  applications  of  the  fluid  dynamics. 

The  kinetic  approach,  is  such  an  alternative,  in  which  we  introduce  some  “smoothed” 
values  of  the  “real”  fluid  parameters,  in  order  to  reduce  the  number  of  degrees  of  freedom  to  just 
a  few  and  to  be  able  to  describe  its  micro-dynamics.  This  assumption  ignores  all  the  system 
characteristics  below  a  certain  threshold  of  the  time  and  space  scales,  and  recover  the  “mean” 
behavior  of  the  fluid,  in  the  statistical  meaning,  after  averages  on  sufficiently  long  time  and  over 
large  enough  microscopic  ensembles. 

The  basic  equations  of  the  kinetic  theory  give  the  evolution  equation  of  f(t,  x,  Q),  the  one- 
particle  phase-space  distribution  function,  in  the  presence  of  the  collisions.  The  distribution 
function  f(t,  x,  £2),  gives  the  complete  statistical  description  of  the  atomic  “ensembles”  and  is 
used  to  define  the  average  values  of  the  fluid.  Usually  the  average  values  are  computed  over 
physical  volumes  dV=  L3,  where  the  characteristic  length,  L,  is  much  larger  than  lm,  the  particle 
mean  free  path,  and  much  smaller  than  Lg,  the  global  space  length  (it  is  quite  obvious  that  we 
have  lm  «  L  «  Lg). 


34 


In  the  absence  of  collisions,  we  have  the  equation  of  conservation  of  the  phase-space,  (the 
Liouville’s  theorem): 

-7-  =  0  (21) 

dt 

where  d/dt  is  a  total  derivative.  For  an  isolated  system  we  can  write  df/dt  as  (using  the  Einstein 
summation  rule): 


f  =d1f  +  vVf=d,f  +  V,d,f  -  02) 

which  represents  the  local  change  in  f  per  unit  time  due  only  to  the  motion  of  particles. 

Now,  if  we  want  to  introduce  the  collision  effects,  the  Liouville’s  equation  becomes: 

^  =  C(/)  (23) 

dt 

where  C(f)  is  a  function  that  models  the  rate  of  change  of  f  caused  by  collisions. 

A  simple  approximation  for  the  collision  operator  was  given  first  by  Boltzmann  as  a  gain- 
minus-loss  operator:  C(f)  =  (G  -L). 

Following  Landau’s  argumentation,  lets  assume  a  two  body  collision  process,  of  particles 
1  and  2,  with  incoming  distribution  functions,  fi0,  f2o>  and  outgoing  distribution  functions,  fic, 
f2c-  If  particle  1  occupied  the  phase-space  dQoi  before  the  collision  and  dQoc  after  collision  (the 
same  for  type-2  particle:  dHo2  is  phase-space  before  the  collision,  and  d£2oc  after  the  collision), 
dQoc  will  not  be  in  dQoi  after  collision  and  the  particle  1  is  said  to  be  lost. 

The  probability  of  loss  will  be  proportional  with  the  number  of  particles  already  in  the 
volume  dV,  namely  fi0,  the  number  of  type-2  particles  that  enter  the  volume  from  the  phase- 
space  range  d£2o2>  namely  g2-df2o2>  the  total  volume  of  allowed  outgoing  phase-space, 
df2icdQ2ci  and  the  probability  of  the  collision  process  Pf  {fi}. 

If  N(t,  x),  is  the  density  function  of  particles  over  all  space,  therefore  [N(t,  x)  dV]  is  the 
mean  number  of  particles  in  the  volume  dV  and  the  total  number  of  losses  L  ,  from  the  phase- 
space  volume  element  dVdQ,  due  to  the  binary  collision  process,  can  be  written: 

L  =  \Pf  {ft}  /,/,  dQ2dQ2cdQlc  dVdQ.  (24) 

Similarly,  particle  gain  into  the  phase-space  volume  d£2  can  only  come  from  the  reversed 
processes,  fic,  f2c t0  flo>  f2o  *  with  fixed  dQqi  and  summed  over  all  d£2ic,  d&2c>  df2o2: 


35 


G  =  J />,<£!)  /u/2t  dQjdQjA 


(25) 


The  Boltzmann  form  of  the  (G  -L)  collision  term  assumes  only  two-body  collisions, 
pairwise  statistically  independent,  with  detailed  (or  at  most  semi-detailed)  balance  symmetry  for 
collision  probabilities. 

From  (23),  (24)  and  (25)  we  finally  obtain  the  Boltzmann  transport  equation: 

(<9,  +  V  id)’f~G~L  -  W 

If  the  system  is  uniform  in  space,  any  distribution  function  f  will  relax  monotonically  to 
the  macroscopic  Maxwell-Boltzmann  form: 


fM_B  =  P  •  exp{-E(p,  v)}  (27) 

where  the  macroscopic  variables  p.v  and  T  (density,  macro-velocity  and  temperature)  are 
independent  of  position  (x). 

In  the  non-equilibrium  case,  any  distribution  function  will  relax  monotonically  in  velocity 
space  to  a  local  Maxwell-Boltzmann  form  and  therefore  p,  v  and  T  will  depend  of  space  and  time. 

In  order  to  recover  the  Euler  equation,  at  the  macroscopic  scale,  from  the  Boltzmann 
transport  equation,  we  use  the  conservation  laws  for  the  collisions  (we  assumed  that  collisions 
preserve  the  conservation  laws  for  mass  and  momentum  exactly,  in  obtaining  the  Boltzmann 
equation): 


\C(f)dQ  =  0 


(28) 


and  respectively: 


J  v  •  C{f)d£l  =  0  (29) 

The  integration  of  the  Boltzmann  equation  will  then  read: 

lU  +  V,a)  fdQ  =  \(G-L)dQ.  (30) 

and  respectively: 

Jv-(a  +  V  rd)-fdn  =  \y-(G~L)dQ.  (31) 
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From  (30)  we  obtain,  after  integration  and  using  (28),  the  continuity  equation: 


[d,P+di(PVi)]  =  ° 


(32) 


where  p(r,x)  =  J/(r,x,£2)dQ  is  the  macroscopic  density  of  the  gas,  and  V  =  iJu-/(r,x,fl>)</£2  is 

P 

the  macroscopic  velocity  vector  of  the  gas,  and  u  its  microscopic  velocity  vector. 

From  (31)  we  will  obtain,  after  integration,  the  momentum  tensor  equation: 

[d,(P\i)+djUij\=o  (33) 

where  the  momentum  flux  tensor  is  given  by: 

nr(v,V;/^  (34) 

As  we  assumed  before,  to  obtain  the  Boltzmann  equation,  in  each  elementary  volume  of 

gas,  dV,  we  have  a  local  Maxwell-B  oltzmann  distribution,  and  in  this  assumption  it  can  be  shown 
from  (34)  that  the  momentum  flux  tensor  has  the  form:  n,>  =pv,V/+5yP>  where  p  is  the 

pressure. 

With  this  expression  for  n,;*  we  obtain  from  (33)  the  same  expression  for  the  Euler’s 
equation  as  in  (12),  from  the  continuum  hypothesis. 

To  recover  the  Navier-Stokes  equation  from  the  kinetic  theory,  as  usual  we  assume  that 
the  gas  reaches  local  equilibrium  in  a  collision  time,  the  one-particle  distribution  function  f(t,  x, 
£2),  has  a  local  Maxwell-B  oltzmann  form,  f|,  and  the  collective  modes  develop  at  large  distances 
and  at  times  much  greater  than  the  molecular  collision  times. 

Following  the  Hilbert’s  approach  (later  used  in  a  some- what  different  form  by  Chapman 
and  Enskog  to  derive  the  transport  coefficients  for  the  macro-dynamic  equations),  it  can  be 
assumed  for  f(t,  x,  Q),  a  perturbation  expansion  of  type  (Hasslacher  1987,  Wolfram  1986): 

f  =  fi(l  +  e(l)  +  em+...)  (35) 

which  can  be  written  also  in  spatial  gradient  expansion  as: 

/  =  /,(!  +  &,  (v)  •  (AV)  +  b2  0 v)  ■ ■  (  XV)2 + . )  (36) 
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where  X  is  the  mean  free  path  (mfp)  of  gas  particles  and  v  is  the  macro  velocity.  The  expansion 
has  to  verify  at  its  n-th  order  the  integral  equation: 

Bn  =  fiC{.eW)  (37) 

where  C  is  the  Boltzmann  collision  operator  from  (23)  and  Bn  is  an  operator  that  depends  only  on 
lower  order  spatial  derivatives.  This  generates  a  recursive  relation  forfifa),  whose  solubility 
conditions  at  order  (n)  are  the  (ml)*  order  hydrodynamic  equations. 

If  we  assume  n=l  (the  first  order  expansion  in  £)  will  have: 

/  =  /((  l+£m)  (38) 

and  from  (26)  we  can  write  (keeping  only  the  first  order  expansion  in  the  collision  operator  and 
putting  £(°)  =  0  in  the  streaming  operator): 

A+v.-a+a.—  /,  =  /|C(£(1))  (39) 

which  is  of  the  form  (37):  B\  =  //  C(£^)). 

From  the  solubility  conditions  (the  Euler  equation  for  p ,  v  and  T,  and  the  ideal  gas 
equation  of  state)  of  (39),  Bi  must  be  orthogonal  to  the  five  of  the  zero  eigenmodes  of  C(£(B)  = 
0  (the  solution  are  1,  v  and  v2). 

In  this  way  we  obtain  a  sequence  of  hydrodynamic  equations,  with  explicit  forms  for  the 
transport  coefficients,  and  from  which  the  order  zero  gives  the  Euler  equation,  order  one  gives  the 
Navier-Stokes  equation  and  order  two  and  higher  give  the  generalized  hydrodynamic  equations 
(they  have  some  validity  only  in  some  special  situations).  Solving  explicitly  for  the  various  £(“), 
provides  a  way  of  evaluating  the  transport  coefficients  (e.g.  viscosity)  of  the  macro 
hydrodynamic  equations. 

2.  2.  The  FHP  models 

The  kinetic  approach,  in  its  essential  features,  was  used  to  devise  the  simplest 
deterministic  local  rules,  made  for  a  collection  of  few-bit,  finite-state  machine,  that  has  the 
Navier-Stokes  equations  as  its  macro-dynamical  description  (we  named  here  the  lattice-gas 
cellular  automaton  fluid). 
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In  constructing  the  lattice  model  of  a  fluid,  one  introduces  major  simplifications  (of 
considerable  computational  convenience),  by  discretizing  space  (point  particles  on  a  lattice),  time 
and  velocity.  Each  node  on  the  lattice  will  be  updated  each  time  step  according  to  a  set  of 
collision  rules,  connecting  the  neighboring  nodes  of  the  lattice,  rules  that  satisfy  conservation 
laws  of  mass  (i.e.  particle  number),  momentum  and  energy  (fig.  1). 

The  first  fully  deterministic  lattice-gas  model  with  discrete  time,  positions  and  velocities 
was  introduced  by  Hardy,  de  Pazzis  and  Pomeau  in  1976,  on  a  square  lattice  (HPP),  but  it  has  had 
only  limited  applications  because  its  hydrodynamic  limit  is  anisotropic  (consequence  of  the 
constraints  imposed  by  the  square  lattice).  The  difficulties  of  the  HPP  model  in  coping  with  the 
full  fluid  dynamics  were  overcome  in  the  model  proposed  by  Frisch,  Hasslacher  and  Pomeau 
(FHP)  in  1986  (Frisch  et  al.  1986),  for  the  2D,  and  that  of  d’Humieres,  Lallemand  and  Frisch  for 
the  3D  case  (d’Humieres  et  al.  1986). 

The  FHP  model  (Frisch  et  al.  1987,  Hasslacher  1987,  Wolfram  1986)  is  constructed  of 
discrete,  identical  particles  which  move  from  site  to  site  on  a  triangular  matrices,  colliding  when 
they  meet,  always  conserving  particle  number  and  momentum.  No  more  than  one  particle  may 
reside  at  a  given  site  and  move  with  a  given  velocity  (the  exclusion  principle)  and  thus  each  node 
of  the  lattice  can  be  described  by  a  six-bit  word  whose  ones  represent  particles  moving  with  the 
velocities  associated  with  their  bit  positions  within  the  word. 

Each  discrete  time  step  of  the  lattice-gas  is  composed  of  two  steps:  first,  each  particle 
hops  to  a  neighboring  site,  in  the  direction  given  by  its  velocity,  and  in  the  second  one,  the 
particle  may  collide,  conserving  mass  and  momentum  (In  Figure  1:  Each  arrow  represents  a 
particle  of  unit  mass  moving  with  unit  speed,  one  lattice  unit  per  time  step,  in  one  of  six  possible 
directions  given  by  the  lattice  links). 

The  precise  collision  rules  are  parameters  of  the  model,  and  could  be  deterministic  or  non- 
deterministic.  In  the  deterministic  case,  for  a  head-on  collision  with  input  particles  on  channels 
(i4+3),  there  are  two  possible  pairs  of  output  channels,  such  as  mass  and  momentum  are 
conserved,  namely  (i+1,  i+4)  (see  figure  2a).  We  can  decide  to  choose  always  only  one  of  these 
channels,  and  consequently  we  have  a  deterministic  model.  The  deterministic  models  are  not 
invariant  under  mirror-symmetry. 

Alternatively,  for  the  non-deterministic  models,  we  can  make  either,  a  random  choice, 
with  equal  probabilities  to  restore  mirror-symmetry,  either  a  pseudo-random  choice,  dependent 
on  a  specific  parameter  (e.g.  parity  of  time  or  space  index,  etc.). 

The  head-on  collision  rules,  conserve  in  addition  of  the  mass  and  momenta,  the  difference 
of  particle  numbers  in  any  pair  of  opposite  directions  (i,  i+3),  which  gives  a  total  of  four  scalar 
quantities.  It  means  that  in  addition  of  mass  and  momentum  conservation,  we  have  a  spurious 
conservation  law.  The  consequence  is  that  the  large-scale  dynamics  of  such  a  model  will  differ 
drastically  from  ordinary  hydrodynamics,  unless  the  spurious  conservation  law  is  removed.  One 
simple  way  is  to  introduce  triple  collisions:  (i,  i+2,  i+4)  — >  (i+1,  i+3,  i+5)  (see  figure  2b). 
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FIGURE  1.  Two  successive  microscopic  configurations  in  the  typical  cellular  automat 
model  discussed  in  Section  2.  Each  arrow  represents  a  discrete  "particle"  on  a  linl 
hexagonal  grid.  Continuum  behavior  is  obtained  from  averages  over  large  numbers  of  pa 


FIGURE  2.  Collision  rules  for  the  FHP  models:  (a)  head-on  collision  with  two  output  channels 
given  equal  weights;  (b)  triple  collision;  (c)  dual  of  head-on  collision  under  particle-hole 
exchange;  (d)  head-on  collision  with  spectator;  (e)  binary  collisions  involving  one  rest  particle 
(represented  by  a  circle). 
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PARTICLE  DIRECTIONS 
IN  THE  HEXAGONAL  MODEL 

FIGURE  3.  The  velocity  vector  of  each  particle  can  point  in  one  of  six  possible  directions.  All 
particles  have  the  same  speed  c. 
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Several  FHP  models  has  been  defined  on  a  triangular  lattice  (Frisch  et  al.  1987, 
McNamara  and  Zanetti  1987,  Chen  et  al.  1992,  Qian  et  al.  1992,  Rothman  and  Keller  1987).  The 
first  one  (FHP-I)  involves  the  simplest  set  of  collision  rules  with  no  spurious  conservation  law: 
pseudo-random  binary  head-on  and  triple  collisions.  FHP-I  is  not  invariant  under  duality 
particle/hole  exchange,  but  can  be  made  so  by  inclusion  of  the  duals  of  the  head-on  collisions 
(see  figure  2c).  It  can  be  further  completed  by  inclusion  of  the  head-on  collisions  with  “spectator” 
(particle  which  remains  unaffected  in  a  collision,  see  figure  4d). 

A  second  type  of  FHP  model  (FHP-II),  is  a  seven-bit  variant  of  the  FHP-I,  including  a 
zero-velocity  “rest  particle”  (see  figure  2e  for  the  additional  collision  rules,  and  figure  3,  for  the 
configurations  of  the  seven  bit  model).  Binary  collisions  involving  rest  particles  remove  spurious 
conservation,  and  do  so  more  efficiently  at  low  densities  than  triple  collisions.  Finally,  the  FHP- 
TTT  model,  is  a  collision-saturated  version  of  the  FHP-II  model . 

The  dynamics  of  the  FHP  models  are  invariant  under  all  discrete  transformations  that 
conserve  the  triangular  lattice:  discrete  translations,  rotations  by  jc/3,  and  mirror-symmetries  with 
respect  to  a  lattice  line. 

There  are  two  main  ways  in  deriving  the  usual  fluid  dynamics  equations  from  the  discrete 
lattice-gas  models  structure:  One  starts  from  the  kinetic  theory  approach  (Wolfram  1986, 
Hasslacher  1987,  Rothman  and  Zaleski  1994)  and  concentrates  on  average  particle  distribution 
functions,  and  the  other,  equivalent  approach,  concentrates  on  the  microscopic  distribution 
functions  (Frisch  et  al  1987,  Rothman  and  Zaleski  1994).  In  the  following  the  main  aspects  of 
both  approaches  will  be  emphasized. 

2. 2. 1.  The  kinetic  equations  of  lattice-gases 

The  lattice  form  of  the  kinetic  single-particle  distribution  function  laws,  can  be  obtained 
from  the  continuous  expressions  by  making  the  following  identifications  between  the  continuous 
and  discrete  probabilistic  formalism: 


/(r,x,Q)->/i(x,r) 

(40) 

|  /(t,x,f2)  dQ.  E,/;(x^)  =  P 

(41) 

Pv-^X,ci-/.(x,f)  =  pu 

(42) 

where  fi  is  the  discrete  single-particle  distribution  function,  giving  the  probability  of  finding  a 
particle  with  velocity  cj  at  position  x  and  moment  t  The  q  are  the  components  of  the  particle 
unit  speed  on  the  directions  given  by  (see  figure  4): 
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FIGURE!  4.  List  of  configurations  for  the  seven-bit  models.  Configurations  with  four  particles 
and  more  are  obtained  by  duality  replacing  particles  by  holes  and  holes  by  particles. 
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,  i  =  1, 2, 


,6. 


(43) 


With  these  notations,  the  lattice-gas  equivalent  of  Liouville’s  theorem  (21)  from  the 
continuum  kinetic  theory  (conservation  of  the  probability  in  the  absence  of  collisions),  can  be 
written  as: 


y.(x  +  h,r  +  k)~  f.(x,t)  =  0 


(44) 


where  h  =  q  dx ,  k  =  dt,  and  dx  ,  dt «  1. 

If  we  expand  the  first  term  of  the  equation  (44)  out  to  02(h,  k)  using  the  Taylor  series 
expansion: 


/,-(xo+h,f£,+*) 


Y'-Uhdx+kStfif  (xo'to)+0(h3,k3) 

x=o  A !  ,• 


(45) 


we  obtain : 


dtdtfi+dxciVfi+^d}dJfj+2t&  (cj- V)2  dt  f  i +dxdt  (c,-  •  V)  £*  /,- = 0 


(46) 


and  to  the  lowest  (first)  order  in  h  and  k  we  have: 


3,f,+  C,-Vf=0  (47) 

which  is  the  lattice-gas  analogous  of  (22),  the  kinetic  transport  equation  in  the  absence  of 
collisions,  which  expresses  the  conservation  of  fluid  (the  continuity  equation).  It  is  the  first 
example  of  macroscopic  equation  for  the  average  behavior  of  a  cellular  automaton  fluid  and  it 
implies,  as  expected,  that  fi  is  unaffected  by  particle  motion  in  a  spatially  uniform  system. 

Now  we  can  write  the  full  Boltzmann  equation,  including  the  collisions,  similar  to  (23) 
from  the  kinetic  theory: 


<?,/,+ <vV/,=C,(/>  <48> 

where  Ci(f)  is  the  discrete  lattice-gas  equivalent  of  the  Boltzmann  collision  operator.  It  gives  the 
time  evolution  of  f,  in  terms  of  two-particles  and  higher  order  distribution  functions  that  appear  in 

Cj(f). 

The  distribution  function  f}  typically  determine  the  macroscopic  average  quantities. 
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In  a  uniform  system  we  have  V  f.  =  0  and  from  (48)  will  have: 

dtfrO  (49) 

Accordingly  we  can  write  the  lattice-gas  form  of  (28)  and  (29)  from  the  kinetic  theory 
(with  the  usual  conventions  (40),  (41)  and  (42))  as: 

Lc,(/)  =  0  (50) 

and  respectively: 

Lc,C(/)  =  0  (51) 

Further  on  we  have  (again  just  following  the  path  from  the  kinetic  theory)  the  continuity 
equation  and  the  momentum  conservation  equation,  that  follow  from  (40)-(41),  (49)  and  (50)  as: 

d(p+da(pua)  =  0  (52) 

which  is  the  usual  continuity  equation  for  the  lattice-gases. 

Momentum  conservation  yields  respectively: 

aLc,/,  +  Lc,(c,-V/,)  =  0  (53) 

where  we  now  define  the  momentum  flux  tensor  (Landau  and  Lifschitz  1959)  as: 

IL*  =  'LiCiaCipfi  (54> 

where  the  Greek  indices  label  the  space  dimensions.  The  momentum  stress  tensor  II  ^  has  to  be 

isotropic  up  to  order  4  in  order  that  the  leading  terms  in  the  momentum  equation  (corresponding 
to  the  convective  and  viscous  terms  in  the  Navier-Stokes  equation)  be  isotropic. 

Finally  the  momentum  conservation  equation  becomes: 

dt(PVct)+  daTla0  =  Q 

For  local  equilibrium  the  macroscopic  distribution  f*(x,  t)  depends  only  on  v(x,  t)  and 
p(x,t),  and  assuming  that  they  vary  slowly  in  space  and  time  (which  is  quite  normal  for 
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hydrodynamic  processes),  and  in  sub-sonic  limit  (||u| «  1),  the  distribution  function  fj(x,  t)  can 
be  approximate  by  a  series  expansion  in  macroscopic  variables  (Chapman-Enskog): 

f.  =  fU  +  A c .  ■  u + B [(c.  •  u)2  - i|u|2  J  +  D J(c .  •  VXc.  •  u) -  j V ■  uj+ . j  (56) 


where  A,  B,  C  are  undetermined  coefficients  and  we  kept  only  the  terms  for  the  standard 
hydrodynamic  phenomena  (the  first  three  terms  account  for  the  change  in  microscopic  particles 
densities  as  a  consequence  of  changes  in  macroscopic  fluid  density,  and  the  fourth  term 
represents  the  first  order  dependence  of  the  particle  densities  on  macroscopic  spatial  variations  in 
the  fluid  density).  We  used  also  the  relations  (41),  (42)  and: 


C,„C.a 
ia  ip 


M.5 

D  uaf) 


(57) 


where  M=6  for  a  hexagonal  lattice  and  D  is  the  number  of  space  dimensions  (2  in  this  case). 

From  (41)  and  (56)  we  obtain  immediately  /  —  while  from  (42)  and  (56),  we  obtain  A 

=  2.  The  coefficients  B  and  C  has  to  be  obtain  only  from  the  explicit  solution  of  (48),  including 
the  collision  terms.  For  uniform  equilibrium  systems  with  u  =  0,  all  fi  are  given  by:  fi  =  f  — 

In  this  case,  the  momentum  flux  tensor  (54)  is  equal  with  the  pressure  tensor,  given  as  in 
standard  kinetic  theory  of  gases  by: 

Paf>='Z‘iCi«Cillf‘iPSc'li  (58) 


where  we  used  also  (57)  for  the  second  equality.  It  gives  the  state  equation  relating  the  scalar 
pressure  to  the  number  of  density  of  the  lattice-gas  fluid:  P-%- 

When  u  0,  the  momentum  flux  tensor  can  be  evaluated,  in  the  Chapman-Enskog 
approximation,  using  the  relations: 


CiaCip 


ciy 


=  0 


(59) 


and  respectively: 


U  Ciccaip  ciyci8  ~  D(D  +  2) 


(S  aS  s  +  8  Sax  +  d  8„  > 

v  a/}  y8  ay  fiS  ad  fiy 


(60) 


Finally  we  obtain: 
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(61) 


n*  +f  BK“f  4m2  V+*  < Cld«u?  4Vu| 

Substituting  in  (55)  and  using  the  momentum  and  continuity  equations,  we  obtain  the 
Navier-Stokes  like  equation  for  the  lattice-gas  fluid: 

(9  (pu)  +  ^  p(u  •  V)u  =  -V(ip +ip|u|2)  -  |CV2  u  +  0(u3)  (62) 

The  standard  form  of  the  Navier-Stokes  equation  for  a  continuum  fluid  in  D  dimensions 
can  be  written  as: 

c?  (p  u)  4-  flp  (u-V)u  =  -Vp+77V2u  +  (£  +  ]|jT7)V(V-u)  (63) 

where  p  is  the  pressure,  T}  and  g  are,  respectively,  shear  and  bulk  viscosity  while  p  from  the 
convective  term,  is  usually  constrained  to  have  value  1  by  Galilean  invariance. 

The  form  of  viscous  term  in  (62)  implies  that  the  bulk  viscosity  g  =  0,  while  the  value  of 
77  is  determined  by  the  coefficient  C  according  to: 

r]=pv=-\pC  (64) 


where  V  is  the  kinematic  viscosity. 

The  convective  term  in  (62)  has  the  same  structure  as  in  the  standard  Navier-Stokes 
equation  but  include  the  coefficient  C/4,  which  can  be  removed  by  a  simple  rescaling  in  velocity: 
u*=  C/4  u,  where  the  coefficients  B  and  C  can  be  obtained  from  the  Boltzmann  approximation. 
From  the  equation  of  state  p  =  p/2  we  obtain  for  sound  speed,  for  the  hexagonal  lattice. 


.  <4  1 

*=7T72 


(65) 


where  D  is  the  number  of  space  dimensions,  and  c  the  lattice  unit  speed. 

For  the  viscosity,  it  can  be  shown  (Wolfram  1986,  Hasslacher  1987,  Rothman  and  Zaleski 
1984)  the  relation: 


V— 


1  1  -1 
12 d(l-d)3  8 


(66) 


47 


where  1/8  is  a  correction  due  to  the  finite  size  of  the  lattice,  and  d=  p/M  is  the  mass  density  per 
cell. 


2. 2.  2.  The  Boolean  kinetic  equations  of  lattice-gases 

As  we  saw,  the  FHP  models  are  defined  on  a  triangular  lattice  (see  figure  3 )  with  unit 

spacing,  and  in  each  node  there  are  six  links  (Ci,  i=l, . 6)  to  its  nearest  neighbors.  Particles  of 

unit  mass  and  velocity  move  along  these  links  in  such  a  way  that  they  reside  in  the  nodes  only  at 
integer  times.  An  exclusion  principle  governs  the  occupation  of  the  links:  no  more  than  one 
particle  can  occupy  any  a  given  time  a  given  link. 

The  state  of  a  node  can  be  describe  as  a  six-bit  Boolean  variable  n  =  {%  i  =  1, — ,  6}, 
where  the  bit  ni  indicates  the  presence  (1)  or  the  absence  (0)  of  a  particle  in  the  i-th  link. 

The  evolution  of  the  system,  from  time  (t)  to  time  (t+1)  can  be  decomposed  in  two  phases: 
collision  and  propagation.  In  the  propagation  phase  each  particle  is  displaced  in  the  direction  of 
its  velocity:  ni(r)  — >  ni(r  +  q).  where  r  denotes  the  sites  of  the  lattice.  The  collision  phase 
conserve  the  mass  and  momentum  (in  order  to  recover  the  Navier-Stokes  equations  on  the 
macroscopic  scale).  In  order  to  avoid  the  presence  of  the  spurious  invariant  (the  difference  of 
particle  numbers  in  any  pair  of  opposite  directions)  triple  collisions  and/or  rest  particles  is 
introduced. 

In  order  to  introduce  the  probabilistic  description  of  the  Boolean  lattice-gas  models,  we 
need  some  basic  assumptions  (Frisch  et  aL  1987): 

i.  Let  s  =  {si,  i  =  1, . .  b}  and  s’  =  {s’i,  i  =  1, . .  b},  denote  the  state  of  a  node  before 

and  after  a  collision  respectively,  with  b  the  number  of  links  per  node,  and  let: 

A  (s  — >  s’)  (67) 

be  the  sit-independent  transition  probability  from  s  to  s’.  The  obvious  normalization  condition  is: 
X,.  A(s  s’ )  =  1,  Vs,  and  the  semi-detailed  balance  condition  then  reads: 

Xs  A(s-»s’)  =  i,  Vs'  (68) 

This  equality  implies  that  the  situation  where  all  states  have  the  same  probability  is  stationary 
with  respect  to  collisions. 

Will  follow  the  approach  from  the  statistical  mechanics:  (i)  definition  of  the  space  phase, 
(ii)  introduction  within  this  set  of  statistical  ensemble  of  initial  conditions  with  a  given 
probability  distribution,  (iii)  evolution  of  probability  distribution  from  the  Liouville  equation,  (iv) 
calculation  of  the  average  values  using  the  evolved  probability  distribution. 
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2. 3.  The  Lattice-gas  Boltzmann  method 


Due  to  the  boolean  dynamics  of  LGA’s  all  simulations  performed  by  this  technique  are 
affected  by  the  statistical  noise,  and  to  get  resonably  resolved  macroscopic  fields  it  is  necessary  to 
avera  ge  over  a  possible  combination  of  large  regions  of  the  lattice,  long  times  and  a  wide  range 
of  initial  conditions. 

In  addition  lattice-gases  usually  are  not  Galilean  invariant  (presence  of  the  so  called  g(p ) 
factor  in  the  nonlinear  advection  term  in  the  Navier-Stokes  equation),  they  have  a  velocity- 
dependent  pressure  and  they  present  spurious  invariants  that  correspond  to  unphysical 
hydrodynamic  quantities. 

The  problem  of  noise  will  disapear  if  in  the  boolean  form  of  the  Boltzmann  equation  (94) 
will  consider  real  variables  Nj  instead  of  the  boolean  ones  (McNamara  and  Zanetti  1988): 

Ni(x+ci,t+\)-Ni(x,t)  =  Ai(N)  (102) 

where  A j(A0  is  obtained  from  the  boolean  collision  term  by  substituting  the  boolean  population 
ni  with  the  ensemble  averaged  population  Nj. 

In  the  double  limit  of  small  Knudsen  and  Mach  numbers  we  can  expand  the  Ihs  of  (102) 
as  (Higuera  and  Jimenz  1989): 


=  A^(p,v)+ATVp,<?v) 

i 


(103) 


and  further  decompose  Neq  as: 

i 


Neq  =  N(0)  +  Njl) + Nf2)  +  0(v3) 


(104) 


where  the  upper  index  refers  to  the  order  in  v,  and  Neq  has  the  form  in  (77). 

i  . 

The  corresponding  expansion  of  the  collision  operator  is: 

A  =  A  ( Nm ) + ^ - iV(1)  +  -^-( N{2)  +  N{neq) ) + _ -4 _ nu) mw 

'■  ^  }  dN.  j  dN{  J  ;  }  2dN  dN  J  7V* 

i 


2dN_dN, 

J  k 


(105) 


where  all  derivatives  are  calculated  at  the  state  of  zero  velocity  Nj  =  d  =  ^ 
As  for  any  equilibrium  distribution,  we  have: 


A.(Neq)  =  0 


(106) 
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and  using  it  for  the  case  of  unifrom  equilibria  (zero  velocity),  one  obtains  also  A f  (1V(0))  =  0 . 
Combining  equations  (105)  and  (106)  one  obtains: 


^M)+A(iv'2) + = o 

dN.  J  dN.  J  >  J  IdNdN,  J  k 

J  J  j  k 


(107) 


Plugging  (107)  into  (105)  we  finally  obtain: 

At(N) -AjjiNj-Np 


(108) 


where  Ajj  is: 


(109) 


or,  introducing  the  transition  probability: 

(no) 

and  p  =  Hisr  The  elements  ay  of  the  collision  matrix  Ay,  are  in  this  case  numerical  parameters 
which  can  be  change  at  will  and  determines  the  scattering  rate  between  directions  i  and  j  (they 
depend  only  on  the  angle  between  the  directions  q  and  Cj). 

The  passage  from  the  complete  collision  operator  to  the  form  (108)  is  clearly 
advantageous  from  the  point  of  view  of  simplicity  and  storage  requirements  of  the  numerical 
scheme. 

The  conservation  of  mass  and  momentum  gives  in  terms  of  the  collision  matrix: 

2XiA,=  o 

(111) 

S?=ic;A/  =  0’  J- . 

For  the  FHP  models  the  number  of  possible  elements  ay  of  the  collision  matrix  is  four  for 
the  6-particle  model  and  we  have  to  add  other  two  for  each  rest  particle  included  (one  to  account 
for  the  influence  of  the  resting  particle  on  itself  and  on  the  other  directions). 
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If  we  denote  by  ae  the  matrix  elements  such  c(-  c.  =  c 2  cos  6  ,  and  as  in  the  hexagonal 
lattice  the  only  possible  angles  are  0,  n/3,  2kJ3  and 7t,  then  (111)  can  be  written  as: 

a0+2a60+2a120+<3180  =0 

(112) 

a0  +  <260-<3120-a180  “° 

As  the  matrix  Ay  is  symmetric  and  cyclic,  the  independent  coefficients  a.Q  can  be 

expressed  in  terms  of  non  zero  eigenvalues  of  the  matrix  (Wolfram  1986),  as: 

X  =  6(a0  +  am) 

(113) 

(7=-6(£Z0+2a60) 


with  multiplicities  2  and  1  respectively.  The  eigenvectors  are  mutually  ortogonal  and  do  not 

depend  on  the  matrix  coefficients,  and  the  eigenvectors  associated  with  the  eigenvalue  X  are  the 
[D(D+2)/2  - 1]  linearly  independent  elements  of  the  set  of  vectors  Qiap- 

To  recover  the  macrodynamics  equations  of  the  in  the  lattice-Boltzmann  method,  one  can 
use  the  same  multi-scale  expansions  around  the  small  parameter  e  (the  inverse  of  the  wavelength 
of  the  typical  spatial  variation  of  the  fields  expressed  in  the  natural  units  of  the  lattice).  Keeping 
only  the  first-order  terms  in  the  multi-scale  expansion,  one  obtains  the  differential  equations: 

dtN{  +(cf  •  d)Ni  =  AfjiNi  -  Nfq)  (1 14) 


In  order  to  better  distinguish  between  the  hydrodynamic  and  non-hydrodynamic  fields,  it 
is  useful  to  project  the  above  eqautions  onto  the  eigenvectors  of  the  collision  matrix  (Frisch  te  al. 

1987).  Due  to  the  symmetries  imposed  by  the  colision  matrix,  all  the  eigenvalues  and  five  of  the 
six  eigenvectors  (1,  cix,  qy  and  the  two  independent  components  of  Qa^)  of  the  collision  matrix 

are  already  known.  One  can  shwon  that  the  other  eigenvector  is  (-1)'.  Therefore,  the  generic 
population  Ni,  can  be  decomposed  as: 

V,-  =  ^P+^pcia  ■  va +| Qaf  Sap +i(-l )'>  (115) 

To  project  equation  (114)  onto  the  orthogonal  basis  of  the  eigenvectors,  one  can  multiply 
successively  by  1,  cia,  Qiap,  (-1)1  and  summing  over  i  will  obtain: 
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(116) 


dtP  +  daJa-0 
d^a  da(P  /  2)  +  dp Sap  —  0 

d^ap  +  Hap  +  RapydyP  =  ^ap~  $ ap  ) 

^t^ap  +4RaPy^ySaP  ~P& 

with: 

Rapy  ~  "5  ^  Qio.pciy 

Saqp=P8(P)(vaVp-^v25ap) 

Hap  =  daJp  +  dpJ a~(dyJy)8ap 

8(p)  =  (3~p)  /  (6-p) 

The  non-isotropy  introduced  by  the  rotational  discretization  of  the  lattice  affects  the 
expression  of  the  non-hydrodynamic  fields,  through  the  ghost  field  ji  and  the  higher-order 
contributions  to  the  stress  tensor. 

In  order  to  recover  the  Navier-Stokes  equations  from  (116),  one  have  to  use  the 
incompressibility  condition  (p~  const.),  the  adiabatic  approximation  (the  time  derivative  term  in 
equation  of  stress  tensor  Sap  be  much  smaller  than  the  other  terms  in  the  same  equation),  and 

finally  to  analyse  the  contribution  of  the  different  fields  in  the  above  equations. 

From  (116)  one  can  observe  that  the  various  fields  have  the  following  order  of  magnitude: 

S = v2  +  ev 

77=£V2+£2V  (118) 

jl  =  O(eri) 

and  as  the  corrections  at  (1 16)  are  at  least  of  second  order,  the  rates  T]/ S  and  fi  IS  tend  to  zero 
as  e— >0  and  thereforethe  ghost  fields  can  be  neglegted  with  respect  to  the  hydrodnamic  fields. 

Substituting  the  stress  tensor  corresponding  to  these  assumptions  in  the  current  equation 
from  (116)  one  obtains  the  Navier-Stokes  equations,  apart  from  the  usual  constant  factor  g(p) 
which  can  be  easily  rescaled  out. 

2.  4.  The  Lattice-gas  Boltzmann  /  Bhatnagar-Gross-Krook  method 

The  lattice-gas  Boltzmann  /  Bhatnagar-Gross-Krook  (LBGK)  method  is  a  further 
abstraction  of  the  initial  concept  of  lattice-gas  models.  As  in  Lattice-gas  Boltzmann  method,  the 
population  densities  are  real  numbers  (and  not  boolean  variables  as  in  the  original  LGA  models), 
the  statistical  noise  is  completely  suppressed  and  the  spurious  invariants  are  easily  controled. 
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However,  what  is  new  in  the  LBGK  method  is  that  the  problem  of  velocity-dependent 
pressure  is  resolved  by  using  a  Maxwellian  distribution  function  that  violates  the  semi-detailed 
balance  of  binary  collisions,  instead  of  the  Fermi-dirac  one  used  in  the  classical  Lattice- 
Boltzmann  method,  and  introduces  more  populated  rest  particles  (Qian  et  al.  1992,  Chen  et  al. 
1992). 

The  Boltzmann  kinetic  equation  was  written  as  in  the  relaxation  method  from  the 
computational  fluid  dynamics  as: 


Ni(x+ci,t+i)  =  (i-T)Ni(x,t)+T  Nie  (x,  t)  -  (119) 

where  Nj  is  the  density  function  of  particle  i,  q  its  velocity  and  Tthe  relaxation  parameter  (from 
the  computaional  fluid  dynamics  the  relaxation  parameter  Thas  to  be  between  0  and  2  to  ensure 
the  stability  of  the  method). 

The  relaxation  proces  is  used  in  order  to  replace  the  collision  term  (without  changing  the 
propagation  term). 

The  equilibrium  distribution  (up  to  the  second  order  to  have  the  correct  nonlinear  term  in 
the  Navier-Stokes  equation)  Nie,  is: 


u  ua  c.  c.a 

_l_  a_J3_/_ta_tP 


2c 


M 


(120) 


where  (X  and  /?  represent  the  Cartesian  coordinates  (with  implied  summation  for  repeated 
indices),  cs  is  the  speed  of  sound,  the  index  p  is  the  square  modulus  of  particle’s  velocity,  and  tp 
is  the  corresponding  equilibrium  distribution  for  u  =  0.  The  tp’s  are  determined  to  achive  the 
isotropy  of  the  fourth-order  tensor  of  velocities  and  Galilean  invariance  (Qian  et  al  1992). 

The  difference  from  the  classical  Lattice-Boltzmann  method  is  that  the  collision  matrix 
Ajj  is  now  replaced  only  by  a  single  value  T. 

Using  the  the  multi-scale  technique  from  the  lattice-gas  models  one  can  obtain  the  Navier- 
Stokes  equations  at  the  second  order  of  approximation: 


d,p+da(pua)  =  0 

^ ci (pua)+ da(puaiip )  d a (pcs  )  +  vdp [ dfj ( pu a )  +  da (p Up )] 


(121) 


1  1  o 

where  cs  —  ■—=  is  the  speed  of  sound,  and  V  =  -i(-|r — l)is  the  viscosity  and  there  are  independent 

of  space  dimension.  For  small  T‘s  the  finite  size  of  the  lattice  gives  a  big  Knudsen  number  and 
the  results  diverges  from  the  the  previous  ones  because  the  system  now  describes  a  rarefied  gas 
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instead  of  hydrodynamics.  It  is  remarkable  that  by  a  proper  choice  of  the  equilibrium  distribution 
the  lattice-gas  BGK  method  was  able  to  simulate  the  fluid  dynamics  (as  it  leads  to  Navier-Stokes 
equations)  and  it  suggests  that  through  a  properly  chosen  equilibrium  distribution  a  wide  range  of 
non-linear  partial  differential  equations  can  be  simulated  in  a  simple  and  efficient  way. 


2. 5.  The  3D  lattice-gas  models 

Three  dimensional  regular  lattices  do  not  have  enough  symmetry  to  ensure  macroscopic 
isotropy  and  therefore  the  way  out  of  this  problem  was  to  move  to  a  higher-dimensional  lattice, 
the  face-centered  hypercubic  (FCHC)  lattice  (figure  5). 

The  nodes  (x^  X2,  X3,  X4)  of  the  lattice  satisfy  the  the  condition:  xi  +  X2  +  X3  +  X4  =  even, 
where  Xj’s  are  integers  numbers.  In  each  node  there  are  24  links  to  its  nearest  neighbours 
(indicated  as  q,  i  =1,  ....,  24)  of  length  V2 .  Propagation  and  collision  are  defined  in  the  same 
spirit  as  in  FHP  models. 

In  the  applications  the  3D  hydrodynamics  periodic  boundary  conditions  are  imposed 
along  the  X4  direction  in  a  layer  of  thikness  2  (along  the  X3  direction)  The  thik  black  line  in  figure 
5  have  a  different  weight  in  the  sense  that  two  particles  can  propagate  along  these  directions, 
according  to  the  two  values  of  the  fourth  component  of  the  velocity  (q)4  =  ±  1. 

It  can  be  shown  (Frisch  et  al.  1987)  that  the  conserved  fourth  component  of  the 
momentum  behaves  as  a  passive  scalar  not  influencing  the  other  three  (there  is  not  a  spurious 
invariant). 

Consider  a  FCHC  lattice  in  a  space  of  D  =  4  dimensions,  with  bm  particles  of  unit  mass 
per  node,  moving  with  velocities  q  (i  =  1,  ...,  bm).  The  component  a  ( a=  1,  ....,  D)  of  these 
velocities  is  denoted  by  cia,  and  the  norm  ||c||,  by  c.  Mc  particles  of  unit  mass  and  zero  velocity 
are  also  present  at  each  node  of  the  FCHC  lattice. 

If  Nj  (r,  t)  (i  =  1,  ...,  bm)  is  the  mean  population  of  moving  particles  at  the  node  r  and 
time  t,  and  N0(r,  t)  is  the  mean  population  of  each  rest  particle,  the  as  usual  we  can  write  the 
Boltzmann  equation  for  the  lattice-gas  as: 

Ni(r+ci,t  +  i)  =  Ni(r,t)  +  Ai[Nl(T,t),N2(r,t), . ,Nb  (r,t),McN0(r,t)]  (122) 

m 

with  i  =  1, ....,  bm,  which  enable  us  to  calculate  the  mean  populations  at  time  step  (t+1),  provided 
that  the  function  A,-,  is  known. 

The  usual  conservation  relations  for  mass  and  momentum  read: 

j  AT  (r + q ,  r + 1)  +  MCN0 (r,  t  + 1)  =  X?=  \  N((r,t)+  MCN0 (r, t )  ( 1 23) 
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v4  =  0 

v4  =  +1  ond  -1 


FIGURE  5.  The  pseudo-four-dimensional  FCHC  model.  Only  the  neighborhood  of  one  node  is 
shown.  Along  the  dotted  links,  connecting  to  next-nearest  neighbors,  at  most  one  particle  can 
propagate,  with  component  V4  =  0;  along  the  thick  black  lines,  connecting  to  nearest  neighbors, 
up  to  two  particles  can  propagate,  with  components  V4  =  ±1. 
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and  respectively: 


X?  =  i  c^-Cr+Cf.f + 1)  =  Xf=  i  C^,(r,r) 


(124) 


The  macroscopic  quantities  p  and  pu ,  are  related  with  these  mean  populations  by: 


p(r,t)  —  2/t  i  A/,(r,f)+McA^0(r,r) 
pU(r,r)  =  2f=  i  CjiVj(r,0 


(125) 


At  equilibrium,  it  can  be  shown  (Frisch  et  al.,  d’Humieres  et  al.  1987)  that  the  mean 
populations  are  given  by  a  Fermi-Dirac  distribution: 

N*q(r,t)  =  [l  +  exp  (h + q  •  c,.)]-1  (126) 

From  (126)  and  (125)  the  mean  population  can  be  expressed  in  terms  of  u  and  p  and  the 
unknown  functions  can  be  expanded  as  a  Taylor  series  around  u  =  0. 

Finally  after  the  standard  Chapman-Enskog  expansion,  rescaling  and  the  multi-scale 
formalism,  the  mean  population  Nj  (r,  t)  can  be  written  as: 


AT.  (r,r)  =  N°  (r,t)+  sNJ  (r,r)  +  £2Nf  (r,  r)+. 


(127) 


where  e=  1/L,  with  1  the  lattice  unit  and  L  the  domain  dimension. 

From  the  conservation  eqautions  and  (127),  up  to  a  second  order  expansion  one  will 
obtain  the  Navier-Stokes  equations: 


<9,p+V-(pu)=o 

5,(p  u) + V  •  P = V  •  [  v(p)V(pu)] + V{[  v(p)  + 1  ]  V  •  (pu)} 


(128) 


where  the  momentum  flux  P  is  expressed  as: 

P = Pcl “ g(P)(“)2 [1 + 0. 5 (D - (^)2 )]  1} + g(p)pu  U  (129) 

S 

The  so-called  galilean  factor  is  given  by: 
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(130) 


g(p)  = 


m\-2d) 

(D+2)bm(l-d) 


and  the  shear  viscosity  v  and  the  bulk  vicosity  £  can  be  written  as: 

b  c 2  c2 

v= - 2 - w - 

D{D+iy  2(Z>+2) 


b  c 2  i  c2  2\ 


(131) 


The  coefficients  y/  and  %  are  derived  from  the  collision  operator  At  once  the  collision 
rules  are  established. 

The  Boltzmann-method  will  read  for  the  FCHC  case  (Gunstensen  and  Rothman  1991): 

Nj (r  +  ci,t+i)  =  N( (r,r)+ X71 "  0  Aj  (r’ r) “  N0j 0] Wy  ( 1 32) 

[  1,  for  j±  0 

with i=l,...,bm, and  W,=  \  . 

1  [Mc,  for  y  =  0 

The  collision  matrix  A  can  be  decomposed  as  follows: 

c°  b°  .  b° 

b° 

:  K) 

b° 


The  matrix  A’y  coprresponds  to  the  collisions  between  the  particles,  while  the 
coefficients  c°  and  b°  correspond  to  collisions  between  moving/rest  particles  and  rest/rest 
particles  repectively. 

Due  to  the  isotropy  of  the  lattice,  the  coeficients  of  A  can  only  depend  on  the  angle 
between  the  velocity  vectors  q  and  Cj  and  on  the  FCHC  lattice  only  5  coefficients  are  possible: 
ao,  a^o,  a.90  and  aiso-  At  these  we  have  to  add  another  two  b°  and  c°  for  the  rest  particles. 

One  can  compute  all  these  coefficients  in  the  same  way  has  been  done  for  the  2D  Lattice- 
Boltzmann  method. 
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3.  Numerical  experiments 


Simulations  and  computations  have  been  performed  with  lattice-gas  models  to  check  the 
basic  properties  of  lattice-gases,  the  linear  response  to  small  perturbations  and  the  behavior  of  the 
transport  coefficients. 

In  the  most  usual  simulations  (e.g.  the  2D  FHP  model  with  rest  particles)  the  state  of  the 
system  at  time  t  is  given  by  a  Li  X  Li  matrix  of  6-7  bit  words  assigned  to  each  node,  and  updated 
regularly  by  collision  +  propagation  steps. 

Boundary  and  initial  conditions  are  set  according  to  the  problem  studied.  For  the  wave 
propagation  experiment,  as  initial  conditions  is  used  a  uniformly  random  distribution  of  particles 
and  velocities,  while  as  boundary  conditions  the  periodic  ones  are  used,  in  order  to  confine  the 
system  on  a  torus  (particles  escaping  at  one  boundary  are  injected  at  the  opposite  one) 

For  the  flow  experiments,  requires  as  initial  conditions,  a  biased  velocity  distribution 
along  a  given  direction,  and  boundary  conditions  that  assure  the  steady  flow  of  incoming  particles 
at  the  input  side  and  constant  density  condition  at  the  output  boundary.  In  the  experiments  with 
solid  boundaries,  depending  on  the  specific  problem,  one  have  to  chose  from  reflection 
conditions  corresponding  to  free-slip  boundaries  (specular  reflection,  see  figure  6a),  no-slip 
boundaries  (bounce-back  reflection,  see  figure  6b)  or  rough  surfaces  (combination  of  specular 
and  bounce-back  reflections  with  equal  probabilities,  see  figure  6c).  Of  course,  the  dimensions  of 
the  obstacles  has  to  be  much  smaller  than  size  L  of  the  space  scale  to  avoid  numerical  artifacts, 
but  that  in  turn  implies  large  number  of  particles  and  therefore  bigger  computer  resources. 

Typical  2D  lattices  are  of  order  of  3  x  106  nodes  (1024  x  3027)  populated  with  6  x  106 
particles,  i.e.  a  density  d  =  0.2.  Stream  line  maps  are  obtained  by  representing  the  velocity  field 
vectors  associated  to  the  fluid  elements  (for  Boolean  LGA’s  by  averaging  the  particles  velocities 
over  a  number  of  nodes:  8  x  8,  32  x  32,  etc.  depending  of  the  problem). 

The  usual  restrictions  regard  the  space  size,  the  minimal  kinematic  viscosity  and  the 
velocity  u  which  must  be  small  compared  to  the  upper  limit  c. 

3. 1.  Poiseuille  Flow 

Poiseuille  flow  as  first  example  of  quantitative  comparison  between  lattice-gas  flow  and 
classical  fluid  dynamics  results,  for  a  system  involving  both  viscous  dissipation  and  non-linear 
behavior  (d’Humieres  andLallemand  1986). 
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a 


b 
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FIGURE  6.  Boundary  reflections:  (a)  free-slip;  (b)  no-slip;  (c)  combination  of  (a)  and  (b). 


FIGURE  7.  Velocity  component  profile  in  a  channel  (a)  close  to  the  inlet,  and  (b)  further 
downstream  [relative  distances  from  inlet  are  ~.5  (a)  and  ~.6  (b)]  (dHumieres  and  Lallemand, 
1986). 
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The  results  presented  here  was  obtained  on  a  2D  FHP  lattice-gas  model  with  512  x  3072 
nodes,  filled  with  a  gas  with  density  d  =  0.22  and  uniform  velocity,  along  the  horizontal  direction, 
u  =  0.3. 

The  system  is  periodic  in  the  vertical  direction  and  wind-tunnel  conditions  are  considered 
at  the  left  and  right  boundary.  The  Reynolds  number  in  the  channel  was  much  smaller  than  that 
corresponding  to  the  appearance  of  turbulence,  and  steady  state  was  attained  after  a  large  enough 
number  of  iterations. 

In  figure  7  the  velocity  profiles  are  presented  for  a  region  close  to  the  input  boundary 
(Fig.  7a)  and  for  a  region  located  about  ten  times  further  downstream,  where  a  characteristic 
Poiseuille  profile  has  developed. 

The  agreement  between  the  lattice-gas  simulations  and  the  profiles  computed  by  the 
Slichting  method  was  good  and  proved  that  the  lattice-gas  models  could  be  an  alternative  at  the 
extremely  laborious  classical  methods  of  fluid  dynamics. 

3. 2.  Flow  past  a  barrier  (von  Karaman  streets) 

Figures  8  and  9  presents  simulation  results  of  a  flow  past  a  plate  (also  among  the  first 
simulated  by  the  lattice-gas  models).  This  2D  flow  is  forced  by  injecting  particles  at  the  left 
boundary  of  the  fluid  space  and  removing  them  at  the  right  boundary  and  thus  creating  a  pressure 
gradient.  The  flow  ,  at  Reynolds  number  of  approximately  300,  creates  vortices,  known  also  as 
von  Karaman  streets,  behind  the  plate,  this  flow  qualitatively  matches  those  that  would  be 
obtained  from  quasi-two-dimensional  experiments  or  other  methods  of  numerical  simulation. 

The  lattice  was  2816  x  1024  nodes  with  periodic  conditions  along  the  y  ax  and  wind- 
tunnel  condition  at  the  left  and  right  edges.  Initially,  the  lattice  is  filled  with  a  gas  of  uniform 
density  and  speed,  with  d  =  0.3  and  u  =  0.428. 

The  presence  of  the  plate  first  produces  shock  waves  due  to  the  reflection  of  the  particles 
at  the  surface  of  the  plate  (figure  8),  then  eddies  start  to  develop  symmetrically  on  either  edges  of 
the  plate.  When  time  reaches  sufficiently  large  values,  it  is  found  that  the  symmetry  of  the  flow  is 
broken  and  vortex  shedding  by  the  plate  occurs,  generating  a  two-dimensional  von  Karaman 
street  (figure  9). 

The  fluid  flow  past  a  cylinder,  for  moderate  Reynolds  numbers  (R  <  100),  was  simulated 
with  a  lattice-Boltzmann  method  for  the  case  of  homogenous  forced  turbulent  flow  (see  figure 
10).  The  results  were  compared  against  a  pseudo-spectral  simulation  for  the  same  case.  Three 
levels  of  resolution  were  adopted:  low  (64  x  64),  moderate  (128  x  128)  and  high  (512  x  512). 
Were  examined  the  range  of  Reynolds  numbers  from  10  to  80  and  were  covered  both,  the  steady 
flows  (Re  <  Rec  =  46)  and  non-steady  periodic  flow.  As  boundary  conditions,  at  the  upstream 
edge  the  mean  populations  for  the  incoming  directions  are  set  to  the  equilibrium  values 
corresponding  to  a  prescribed  value  of  the  uniform  free-stream  velocity.  The  free  stream  velocity 
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FIGURE  8.  Map  of  the  flux  of  particles  in  a  886.8  x  2048  channel  3000  time  steps  after  the 
introduction  of  a  flat  plate  of  real  size  216.5.  The  distance  between  the  back  of  the  plate  and  the 
point  where  jx  =  0  on  the  axis  of  the  channel  is  defined  as  the  size  of  the  wake. 
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FIGURE  9.  Similar  to  Figure  8,  but  after  40000  iterations  showing  the  formation  of  a  Karman 
street.  To  emphasize  the  vortices,  the  mean  momentum  has  been  subtracted  from  the  local  ones. 


was  0.2  (Mach  number  =  0.3)  to  control  the  compressibility  effects  and  spurious  terms  in  the 
momentum  conservation  equation.  At  the  downstream  edge  the  outgoing  populations  are 
multiplied  by  a  factor  smaller  that  1  and  bounced  back,  to  simulate  the  effect  of  a  porous  wall 
with  pores  open  to  vacuum  ( the  reduction  factor  assures  a  uniform  flow  equal  to  the  one  injected 
upstream).  On  the  lateral  edges  of  the  domain,  were  applied  periodic  conditions.  On  the  points  on 
the  polygonal  approximating  the  surface  of  the  body,  were  used  the  bounce-back  conditions. 

In  the  first  stage  of  simulations  acoustic  waves  were  developed  but  after  a  transient  period 
of  8-10  residence  time  D/U  (D  is  the  cylinder  diameter  and  U  is  the  free-stream  velocity)  the 
recirculation  regions  attained  the  steady  state  (Figure  10). 

The  agreement  with  other  numerical  methods  and  experimental  values  was  good  in  the 
range  of  Re  *  40,  and  the  discrepancies  found  for  high  Reynolds  numbers  are  due  to  the  moderate 
size  of  the  computational  domain. 

The  same  problem  has  been  treated  in  three  dimensions  by  Frisch  et  al.  (1988)  using  a  4- 
D  FCHC  lattice  (see  figure  11).  The  initial  and  boundary  conditions  are  roughly  a  3D 
transposition  of  those  used  in  the  2D  simulations  of  flow  past  a  plate.  The  lattice  size  was  128  x 
128  x  256,  with  periodic  conditions  in  the  directions  transverse  to  the  upstream  flow. 

Figure  1 1  gives  an  image  of  the  3D  perspective  view  of  high-vorticity  modulus  regions, 
displaying  a  “basket  and  handle”  structure. 

3. 3.  Channel  flow  in  expanded  geometry 

Another  well-known  flow  situation  is  that  of  two-dimensional  backward  facing  step  at 
low  Reynolds  number  (figure  12).  Here  was  used  a  lattice  of  size  4608  x  512  with  an  inlet  part  of 
length  512  and  width  256  (Noullez  et  al.  1986,  d’Humieres  and  Lallemand  1990). 

The  lateral  boundaries  of  the  channel  and  of  the  backward  facing  step  are  set  with  the 
stick  condition.  At  the  inlet  was  used  the  wind-tunnel  condition  with  injection  of  particles 
distributed  with  uniform  density  and  a  parabolic  velocity  profile.  The  experiments  were  done  at 
50,  100  and  150  Reynolds  numbers  (the  adjustment  of  the  Reynolds  number  is  done  either  by 
changing  the  velocity  or  the  density,  due  to  the  dependence  of  g(p)  vs.  p). 

After  the  system  reaches  the  steady  state  a  recirculation  zone  was  observed  behind  the 
step  (as  was  observed  experimentally)  and  the  location  of  the  reattachment  point  was  computed 
(see  figure  12b).  The  density  of  gas  was  computed  and  found  to  vary  by  less  than  3%  over  the 
entire  lattice,  indicating  that  the  flow  is  essentially  incompressible. 

As  can  be  seen  from  figure  10,  a  good  agreement  with  the  behavior  of  real  fluids  was 
obtained,  for  different  Reynolds  numbers. 
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FIGURE  10.  Streaklines  of  the  far-field  at  Re  =  52.8  two  shedding  cycles  after  the  transient. 
The  elliptical  shape  of  the  cylinder  is  due  to  different  scales  along  the  horizontal  and  vertical 
axes. 
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t  =  4000  =  6.2  tr 

FIGURE  1 1 .  Perspective  view  of  the  high  vorticity  regions  developed  at  Re  =  190  in  a  3-D  flow 
passed  a  circular  plate  (shown  in  black)  (Rivet,  Henon,  Frisch,  and  d'Humieres,  1988). 
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FIGURE  12.  Channel  flow  with  sudden  expansion.  Isomach  curves  for  velocity  component 
along  the  average  input  flow  direction  at  Re  =  50(a)  and  Re  =  150  (b,c);  a u  =  10%.  Isodensity 
curves  (at  p  —  .95  pmax)  at  Re  —  150  (d).  (Note  2x  expanded  spatial  scale  in  c  and  d.)  Arrow  (in  b) 
indicates  reattachment  point  (Noullez,  Lallemand,  and  d'Humieres,  1986). 
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3.4.  Free  boundaries  flow 


The  lattice-gas  models  proved  to  be  extremely  appropriate  to  study  mixture  of  fluids 
(miscible  or  imiscible)  and  an  example  in  that  respect  is  the  simulation  of  the  Rayleigh-Taylor 
instability  (see  Figure  13,  a  heavy  fluid  penetrates  a  lighter  fluid  layer)  and  Kelvin-Helmholtz 
instability  (see  Figure  14,  two  fluid  layers  moving  in  opposite  directions  with  respect  to  each 
other,  develop,  by  shear  constraint,  a  roll  up  at  the  interface). 

The  Rayleigh-Taylor  instability  was  simulated  on  a  512  x  512  lattice,  initially  set  up 
without  overall  motion,  between  two  rigid  horizontal  plates. 

At  the  initial  time,  a  gravity  potential  is  applied  in  a  such  a  way  that  the  initial  state 
become  unstable,  and  depending  on  the  initial  form  of  the  interface,  an  instability  develops  in 
time. 

If  the  interface  is  planar,  a  distortion  will  evolve,  which  shape  depends  on  the  initial 
microscopic  conditions.  If  the  interface  is  sinusoidal,  the  distortion  amplitude  increases 
exponentially  (periodic  boundary  conditions  is  assumed  in  the  horizontal  direction)  and  after 
some  time  they  are  no  longer  harmonic,  and  the  typical  shapes  are  observed  (figure  13). 

The  Kelvin-Helmholtz  instability  (Clavin  et  al.  1990),  was  simulated  between  two  parallel 
plates  of  length  1024,  separated  by  256  lattice  sites,  set  as  channel  boundaries  (assumed  to  be 
periodic  in  the  direction  of  the  plates). 

Stick  conditions  are  implemented  at  the  boundaries  and  the  initial  conditions  set  a  half  of 
the  channel  with  A  particles  with  +u  velocity  (parallel  to  the  plates)  and  the  other  half  with  B 
particles  with  -u  velocity  (see  Figure  14).  In  time,  an  instability  is  developed  and  A  particle  are 
advected  by  the  flow  into  the  B  zone. 


66 


FIGURE  13.  2-D  lattice  gas  simulation  of  the  Rayleigh-Taylor  instability.  Maps  of  A-particle 
flux  (a)  and  of  B-particle  flux  (b)  after  t  =  1600  (Clavin,  d'Humieres,  Lallemand,  and  Pomeau, 
1986). 
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FIGURE  14.  2-D  lattice  gas  simulation  of  the  Kelvin-Helmholtz  instability.  Average  velocities 
at6  to  right)  and  -u  (right  to  left)  in  lower  half  and  in  upper  half  of  channel  respectively. 

Map  shows  flux  of  particles  of  one  species  (d'Humieres,  Lallemand,  and  Searby,  1987). 
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4.  Applications 


4. 1.  Flow  through  porous  media 

The  lattice-gas  models  are  particularly  efficient  in  complex  geometry’s  as  it  is  the  case  of 
the  flow  through  porous  media.  The  lattice-gas  models,  due  to  their  flexibility  in  dealing  with 
complex  boundary  conditions,  and  due  to  their  intrinsic  particle  tracking  techniques,  can  easily 
manage  the  intricacies  associated  with  flow  in  complex  geometries. 

The  first  lattice-gas  simulations  of  flows  in  porous  media  were  performed  with  2D  FHP 
model  (Rothman  and  Keller  1988)  and  this  pioneering  work  was  subsequently  extended  to  3D, 
using  either  the  FCHC  (Rivet  et  al.  1991)  either  the  lattice-gas  Boltzmann  method  (Gunstensen 
and  Rothman  1991,  Benzi  et  al.  1992). 

The  porous  medium  (either  in  2D  or  3D)  is  modeled  as  a  random  sequence  of  elementary 
blocks,  of  at  least  four  lattice  units  wide.  The  distribution  is  such  that  no  pore  with  cross  section 
narrower  than  a  4  x  4  lattice  units  could  be  found  (the  minimal  size  of  pore  could  be  find  by 
running  preliminary  tests  to  ascertain  under  which  parameter  conditions  -  Re,  Mach  number  and 
pore  size  h  -  the  model  reproduces  the  Poiseuille  flow  in  a  square-section  of  the  pore  channel). 

The  pore  size  has  to  be  as  small  as  possible  in  order  to  have  a  large  granularity  of  the 
porous  medium,  and  in  the  same  time  it  cannot  be  too  small,  because  otherwise  gradients  of  the 
order  of  the  pore  size  develop  and  would  impair  the  convergence  of  the  lattice-Boltzmann 
equation  to  the  Navier-Stokes  equations  (see  figure  15). 

Here  again  the  possibility  to  “tune”  the  particle  mean  free  path  turns  to  be  very  important 
in  order  to  improve  the  convergence  to  hydrodynamics  at  small  scales. 

In  general,  the  emphasis  is  on  the  empirical  investigation  of  the  dependence  of  bulk 
properties  of  flow  on  aspects  of  microscopic  disorder:  asses  the  parameter  range  in  which  Darcy’s 
law  is  fulfilled,  and  measure  the  permeability  of  the  medium. 

4. 1. 1.  Immiscible  multi-phase  flow  through  porous  media 

There  are  two  basic  approaches  to  the  simulation  of  multiphase  flow  through  porous 
media.  The  first  is  to  make  some  simplifying  assumptions  about  the  flow,  such  a  linear  force-flux 
law  (e.g.  Darcy  type),  and  also  about  the  geometric  properties  of  the  porous  medium  (e.g.  a 
network  of  pores  and  throats).  The  second  simulates  the  hydrodynamics  of  the  immiscible  fluids 
through  a  microscopic  model  of  the  porous  medium,  without  any  restrictive  assumptions  about 
the  flow  (Rothman  and  Keller  1988,  Gunstensen  and  Rothman  1992). 

The  lattice-gas  model  of  immiscible  fluids  (ELG)  is  based  on  the  Lattice-Boltzmann 
method  in  which  a  new  parameter  was  introduced  to  account  for  the  type  of  fluid:  “the  color’. 
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FIGURE  15.  Lattice-gas  simulation  of  flow  through  a  two-dimensional  porous  medium 
(Rothman,  1988).  The  fluid  is  forced  from  left  to  right. 
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The  particles  are  colored  either  red  or  blue  (or  other  color  if  there  are  more  than  two  phases)  and 
the  usual  collision  operator  is  altered  to  obtain  the  correct  interfacial  dynamics  between  phases. 

No-slip  boundary  conditions  at  the  solid-fluid  boundaries  and  the  appropriate  surface 
tension  boundary  conditions  are  imposed  at  the  fluid-fluid  interfaces. 

The  particle  configuration  at  a  site  is  described  by  the  variables  r  =  {rj,  i  =  1, 6}  for 
the  red  particles,  and  b  =  {bj,  i  =  1, ....,  6}  for  the  blue  particles.  The  evolution  of  the  particles  is 
as  usual:  first  particles  are  moving  to  the  nearest-neighbor  sites  with  velocity  q,  but  in  the 
collision  step,  it  is  conserved  not  only  the  total  number  and  momentum  of  particles,  but  also  the 
“color”! 

To  achieve  the  surface  tension,  the  configuration  resulting  from  a  collision  depends  on  the 
configurations  at  the  nearest-neighbor  sites  of  the  lattice.  To  simulate  attraction  between  the 
particles  of  like  color,  those  collisions  that  send  red  particles  in  the  direction  where  the  neighbors 
contain  a  relative  majority  of  red  particles  are  favored,  while  blue  particles  are  sent  toward 
neighboring  blue  particles. 

Two  new  vector  fields  are  defined:  the  color  flux  and  the  color  field,  respectively: 


qir(x)Mx)}  =  X,  c,[rt(x)  -  (x)] 

/(*)  =  X,  cfX,  c,  [r;-  (x + c, )  -  bj  (x + c, )] 


(133) 


The  particle  configuration  that  results  from  collision  is  then  that  configuration  which 
maximizes  the  product  (f  •  q  )  such  that  the  number  of  reds,  the  number  of  blues  and  total 
momentum  are  conserved  (Rothman  and  Zaleski  1989).  By  these  rules,  the  adherence  of  the 
original  lattice  (un-colored  LB  method)  to  the  Navier-Stokes  equations  applies  equally  well  to 
single-color  regions  of  ILG  (Rothman  and  Keller  1988,  Succi  et  al  1989). 

The  porous  medium  used  in  the  2D  simulations  of  immiscible  fluids  was  constructed  by 
placing  solid  squares  of  random  size  on  the  lattice.  The  length  of  the  square’s  side  is  uniformly 
distributed  between  16  and  8  lattice  units  when  they  are  placed  on  the  128  x  128  lattice  used  in 
simulations. 

Figure  16  illustrates  the  case  in  which  a  nonwetting  fluid  (black)  is  injected  into  a  porous 
medium  filled  with  a  wetting  fluid  (white).  The  solid  sites  are  gray,  and  there  are  an  average  a 
number  of  4.9  particles  per  site  (both  fluids  have  the  same  viscosity). 

When  the  capillary  number  (Ca  =  [iu/  o,  the  ratio  of  viscous  forces  to  surface  tension) 
Ca«  1  ,  capillarity  dominates  the  flow,  causing  the  invading  nonwetting  fluid  to  follow  the  path  in 
which  the  critical  capillary  pressure  is  the  smallest  (in  figure  16,  Ca  =  10*2  small  enough  that  the 
nonwetting  fluid  is  unable  to  penetrate  the  narrowest  channels). 

Figure  17,  illustrates  the  opposite  case,  in  which  the  wetting  fluid  (black)  invades  the 
same  porous  medium  now  initially  occupied  by  the  nonwetting  fluid.  In  the  same  conditions  as  in 
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Time  r  is  given  in  time  steps. 
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the  previous  experiment,  the  wetting  fluid  fills  the  porous  medium  relatively  quickly,  and 
comparatively  smoothly  without  appearance  of  the  capillary  fingers  as  in  Figure  16. 

Some  other  coloring  schemes  were  proposed,  for  the  spinodal  decomposition  in  porous 
media,  based  on  the  LB  method  (Cyeplak  et  al.  1992,  d’Ortona  et  al.  1994,  Rybka  et  al.  1995). 

This  approach,  use  the  “coloring  scheme”  way  of  the  ILG  but  with  a  drastic  change  of  the 
re-coloring  algorithm  (Cyeplak  et  al.  1992,  d’Ortona  et  al.  1994,  Rybka  et  al.  1995).  As  the  phase 
separation  in  porous  media,  due  to  the  spinodal  decomposition,  depends  crucially  on  the  way  the 
surface  tension  effects  are  implemented,  the  new  re-coloring  scheme  allows  for  a  more  flexible 
control  of  the  interfaces  in  the  kinetic  model. 

The  surface  tension  is  controlled  by  two  microscopic  parameters,  one  determining  the 
magnitude  of  the  jump  in  the  density  and  the  other  the  thickness  of  the  interface.  More,  the  LB 
model  is  not  liniarized  around  the  equilibrium  and  is  Galilean-invariant. 

The  recoloring  scheme  is  defined  by  the  equations: 


rL(x)  = 

b'i(x)  = 


% 

Rt  +  B, 


fi+Pi 
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(Rr+B,)2 
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■fi-h 


Mi 

(Rt  +  B,)2 


COS  (pi 


with : 


fiW  =  fi+Pl\g\cos(2(pi) 


(134) 


where  fi  refers  to  the  initial  state  of  density,  fi  is  the  density  after  redistribution,  <pt  is  the  angle 

with  the  i-th  direction  of  the  lattice  links,  Rti  Bt  refers  to  the  total  red  and  blue  densities  at  a  node, 
and  g  is  the  gradient  vector  at  each  node.  The  parameters  /3j  and  j52  control  the  values  of  the 

surface  tension  and  interfacial  width  respectively.  They  allow  to  avoid  having  an  abrupt  interface 
whose  properties  are  sensitive  to  its  orientation  with  respect  to  the  underlying  lattice. 

The  simulation  was  carried  out  on  a  256  x  256  lattice,  the  fluid  density  was  1/3  and  the 
relative  content  of  red  particle  was  0.5,  with  Re  between  3  and  7. 

Figure  18  displays  the  spinodal  decomposition  for  wetting  fluid  (a)  and  a  neutral  fluid  (b). 
For  longer  times  one  can  see  easily  the  difference  in  the  shape  of  separated  phases.  In  the  wetting 
case  the  red  fluid  sticks  to  the  solid  wall  and  the  domains  show  much  greater  connectivity.  In  the 
neutral  case  the  domains  are  more  spherical  and  less  connected  together. 


4. 1.  2.  Phase  separation 

There  are  many  difficulties  in  numerically  simulating  multi-phase  flow  or  immiscible 
fluids  through  porous  media,  with  conventional  methods.  Usually,  the  multi-phase  flow  through 
porous  media  is  a  mixture  of  water/oil/gas  or  any  other  combination  of  two  of  them.  Phase 
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modelled  by  random  threads,  (a)  wetting  and  (b)  neutral  "solid”  nodes. 


separation,  due  to  spinodal  decomposition,  is  associated  with  first  order  transitions  such  as  liquid- 
gas  transition,  fluid  separation  in  binary  fluid  systems  and  in  binary  metallic  solutions. 

One  approach  starts  from  the  LB-ILG  (Rothman  and  Keller  1988,  Gunstensen  et  al.  1991) 
but  introducing  S  total  components  (  on  2D  hexagonal  lattice  or  4D  FCHC): 


Nf  (x + c„  t + l) - N?(x,  t)=  Af  (x,  t) 
G  =  a  =  0, . b 


(135) 


where  N?(x,t)  is  the  single-particle  distribution  function  for  the  cr  component  and  Af(x,f)  is 
the  collision  term.  The  set  of  vectors  { q,  i=  1, ...,  b}  are  the  velocities  of  the  particles  moving  on 
the  lattice  links  to  the  one  of  the  b  nearest-neighboring  sites,  each  time  step  (c0  =  0  is  associated 
with  the  rest  particles). 

The  collision  parameters  has  the  LB-BGK  form  with  the  corresponding  modification  due 
to  the  S  components: 

Af(x,t)  =  -^-[N?(x,t)-N?('eq\x,t)]  (136) 

cr 

where  Ta  is  the  mean  collision  path  (Qian  et  al.  1992,  Chen  et  al.  1992)  for  the  c-th  component 
and  determines  its  fluid  viscosity. 

N°(eq)(x,t)  is  the  equilibrium  distribution  at  site  x  and  time  t,  with  its  explicit  form 
chosen  as  is  specific  for  the  LB-BGK  models. 

In  order  to  simulate  multiple  components  fluids  or  fluids  with  non  ideal  gas  equation  of 
state,  the  nonlocal  interactions  among  particles  are  incorporated  through  an  interaction  potential: 

V(x,x')  =  G(ja,(x,X)  xP<T(x)vF°'(x' )  (137) 

where  (x,x')  is  a  Green’s  function  controlling  the  sign  (attractive/repulsive)  and  the 
strength  of  interaction,  while  ^(x)  plays  a  role  as  the  effective  number  density  for  the 
component  <r. 

With  this  interaction  potential,  the  LB-BGK  model  exhibits  thermodynamic  phase 
transitions  shown  in  Figure  19. 

The  numerical  simulation  was  done  on  a  2D  256  x  256  hexagonal  lattice,  with  relative 
density  d  =  0.5.  The  lattice  was  initialized  with  a  homogenous  density  distribution  with  only  1% 
random  perturbation. 

Figure  19  displays,  in  gray  scale,  the  evolution  of  the  mass  density  distribution.  As  G  gets 
below  the  critical  value,  the  system  changes  from  a  single-phase  to  a  two-phase  fluid. 
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FIGURE  19.  Phase  transition  in  a  LB  simulation  on  a  256  x  256  lattice  with  a  single  component. 
Shown  are  the  time  evolution  of  the  number  density  distribution  at  t  =  20  (a),  t  =  200  (b),  t  = 
2000  (c),  and  t  =  20000  (d).  The  variation  of  the  density  is  shown  in  gray  scale  with  the 
minimum  in  black  and  the  maximum  in  white.  Since  the  lattice  is  hexagonal,  the  graphics  are 
distorted  by  a  factor  of  V3/2  in  one  direction.  G  =  -0.45,  po  =  1,  and  (p)  =  0.693  =  In  2. 
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It  can  be  seen  the  separation  of  the  different  thermodynamics  phases  and  the  existence  of 
the  surface  tension. 

4. 1. 3.  Viscous  fingering 

Viscous  fingering  is  a  complicate  phenomenon  which  arises  as  a  result  of  instabilities  that 
occur  when  a  dense  viscous  fluid,  filling  a  porous  material,  is  forced  to  be  displaced  by  the  local 
injection  of  a  less  viscous  fluid. 

This  instability  is  analogous  to  the  Saffman-Taylor  instability  observed  in  Hele-Shaw 
cells.  In  porous  media  the  underlying  mechanism  is  basically  two  dimensional  and  so  the 
phenomenon  is  well  suited  for  lattice-gas  investigation. 

The  porous  medium  is  simulated  by  a  random  spatial  distribution  of  static  scattering 
nodes  (3%  of  the  total  number  of  nodes)  on  a  320  x  512  hexagonal  lattice.  The  two  fluids  are 
miscible  and  have  a  viscosity  ratio  of  33. 

The  viscosity  difference  is  set  by  selecting  specific  collision  rules  for  each  fluid:  for  the 
more  viscous  fluid  only  12  efficient  collisions  are  used  among  the  64,  which  in  turn  are  all 
present  for  the  low  viscosity  fluid. 

Particles  colliding  with  the  scattering  nodes  of  the  porous  material  undergo  bounce-back 
collisions,  yielding  momentum  dissipation  at  these  nodes.  During  collisions  involving  particles  of 
both  fluids,  the  particles  are  redistributed  so  as  to  minimize  mutual  diffusion,  such  that  the 
mixing  be  a  slow  process  on  the  hydrodynamic  scale.  Minimizing  the  diffusion  coefficient  D  one 
maximizes  the  Peclet  number  (Pe  =  Ud/D,  where  U  is  the  average  Darcy’s  velocity,  d  is  the  pore 
size,  and  D  is  the  diffusion  coefficient,  is  the  ration  between  rate  of  transport  by  convection  to  the 
rate  of  transport  by  molecular  diffusion),  which  favors  the  development  of  large  body  waves 
modes. 

The  lattice  was  initialized  with  a  slab  of  low  viscosity  fluid  along  the  left  boundary  and 
with  a  highly  viscous  fluid  in  the  rest.  At  initial  time,  both  fluids  are  at  rest  and  have  the  same 
density.  The  right  edge  of  the  lattice  is  continuously  maintained  at  this  constant  density  in  order 
to  simulate  an  infinitely  extended  medium.  Periodic  boundary  conditions  are  imposed  on  the 
upper  and  lower  boundaries.  A  pressure  gradient  is  set  up  and  maintained  by  a  density  gradient 
oriented  to  the  left,  which  pushes  the  less  viscous  fluid  from  left  to  right  through  the  medium. 

Viscous  penetration  then  occurs  and  develops  progressively  into  viscous  fingering  (Figure 
20  and  Figure  21)  with  a  characteristic  wavelength,  the  rate  of  growth  being  proportional  to  the 
porosity  of  the  medium. 

The  experiments  have  shown  that  viscous  fingering  occurs  only  if  there  is  a  sufficiently 
large  viscosity  difference  between  the  two  fluids  and  when  there  is  a  density  gradient  on  the  less 
viscous  fluid,  otherwise  the  system  will  exhibit  only  mutual  or  forced  diffusion. 
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FIGURE  20.  Viscous  fingering  in  2-D  porous  lattice.  Time  is  given  in  lattice  gas  time  step 
units.  Notice  tip-splitting  at  t  =  10000. 
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FIGURE  21.  Simulation  of  viscous  fingering:  a  high-mobility  fluid  (entering  from  below) 
displacing  a  much  less  mobile  fluid  (leaving  from  above).  The  simulation  is  done  on  a  256  x  256 
lattice. 


4. 2.  Transport  of  particle  suspensions 


Simulation  of  transport  of  particle  suspensions  by  fluids  is  of  considerable  interest  for 
industrial  processes,  biological  liquids,  in  oil  reservoir  exploitation,  etc.  The  theoretical 
description  of  these  systems  at  particle  volume  fractions  exceeding  the  dilute  limit  is  considerable 
complicated  by  the  presence  of  indirect  (hydrodynamic)  interactions  between  the  particles.  These 
interactions  result  from  the  velocity  fields  set  up  in  the  suspending  liquid  by  the  relative  motion 
of  solid  particles  (Landau  and  Lifsitz  1967).  The  treatment  of  these  forces  is  complicated  by  the 
fact  that  they  are  of  many-body  and  long-ranged  nature  and  therefore  the 'numerical  simulations 
become  extremely  important  of  studying  the  dynamical  properties  of  particle  transport.  Most 
simulations  algorithms  (Brownian  dynamics,  Stokesian  dynamics,  the  multipole  method)  are 
based  on  the  clear  time-scale  separation  that  exists  between  the  dynamics  of  fluid  and  that  of 
solid  particles.  This  separations  implies  that  the  development  of  the  hydrodynamic  interactions  is 
instantaneous  and  therefore  they  depend  on  the  positions  and  velocities  of  all  particles.  For  this 
reason,  all  these  algorithms  scale  as  the  square  root  or  cube  of  the  number  of  particles. 

In  the  Lattice-Boltzmann/BGK  method  (Behrend  1995,  Ladd  1990, 1991,  1993, 1994a,  b) 
the  state  of  fluid  is  updated  on  a  regular  lattice  while  the  solid  particles  move  continuously  in 
space  and  interact  with  the  fluid  at  a  set  of  special  lattice  nodes.  This  technique  takes  advantage 
of  the  fact  that  the  hydrodynamic  interactions  are  time  dependent  and  develop  from  purely  local 
interactions  at  the  solid-liquid  interface.  Therefore  it  is  not  necessary  to  consider  the  global 
system,  but  one  can  update  one  particle  at  a  time.  The  methods  scale  linearly  with  the  number  of 
solid  particles  and  therefore  allows  far  larger  and  more  significant  simulations  than  those  possible 
with  conventional  methods. 

The  hydrodynamic  interactions  between  solid  particles  are  fully  accounted  for  in  the 
lattice-Boltzmann  methods,  both  at  zero  and  finite  Reynolds  numbers  (Ladd  1994  a,  b). 

The  crucial  part  of  the  algorithm  is  the  mechanism  of  interaction  between  the  solid 
particles  and  the  fluid,  the  so-called  “boundary  rules”,  implemented  at  a  set  of  lattice-sites,  the 
“boundary  nodes”.  The  analytical  developments  of  the  effects  of  the  boundary  rules  are  possible 
only  for  simple  geometries  (such  as  planar  Couette  flow)  while  for  more  complex  boundaries, 
such  in  particulate  suspensions,  the  only  judge  of  the  quality  of  a  boundary  method  is  the 
comparison  of  the  results  of  numerical  computation  with  independent  calculations. 

The  solid  particles  are  defined  by  a  boundary  surface,  of  any  size  and  shape,  which  cuts 
some  of  the  links  between  the  lattice  nodes  (see  figure  22)  and  defines  a  set  of  boundary  nodes 
whose  positions  are  rt,.  At  each  update  of  the  lattice  a  special  rule  at  the  boundary  nodes  is 
implemented  on  the  distribution  functions  n;.  This  boundary  mle  exchange  momentum  between 
the  fluid  and  the  particle  (the  total  momentum  being  conserved)  and  enforces  a  stick  boundary 
conditions  on  the  fluid.  Therefore  the  fluid  velocity  at  the  boundary  nodes  is  matched  to  the  local 
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FIGURE  22.  Geometry  for  quasi-periodic  simulations.  The  figure  illustrates  a  system  that 
closely  approximates  a  periodic  one  under  an  external  flow.  The  lattice  nodes  are  shown  by  solid 
circles  and  the  boundary  nodes  by  solid  squares.  The  arrows  indicate  velocity  directions  c,-  and 
d'  at  the  boundary  nodes.  Periodic  boundary  conditions  are  applied  across  the  planes  indicated 
by  solid  lines;  the  configuration  of  solid  particles  within  the  quasi-periodic  unit  cells  (bounded 
by  dashed  lines)  are  identical.  Macroscopic  flows  can  be  set  up  by  the  planes  of  boundary  nodes 
at  either  end  of  the  system.  With  this  geometry  we  can  set  up  uniform  flow  perpendicular  to  the 
boundary  walls  or  an  approximately  linear  shear  flow  parallel  to  the  boundary  walls.  The 
properties  of  the  central  cell  are  close  to  those  of  a  truly  periodic  system;  it  is  not  necessary, 
apparently,  to  include  more  cells,  although  this  can  be  done  if  required. 
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solid-body  velocity  Ub  •  Ub  and  is  determined  by  the  solid-particle  velocity  U,  its  angular  velocity 
Q,  and  the  position  of  its  center  of  mass  R  (Behrend  1995): 

ub=U+£2x(rb-R)  (138) 

The  momentum  density  exchange  causes  a  local  force  density  to  be  exerted  on  the  particle  at 
node  i*b,  f(  i"b).  The  total  force  and  torque  on  the  particle  are  obtained  by  summing  f(rb)  and 
(rb-R)  X  f(rb)  over  all  boundary  nodes  associated  with  the  particle.  These  forces  and  torque  are 
used  to  update  the  position  and  velocities  of  the  particle  according  to  the  laws  of  Newtonian 
mechanics,  using  pre-assigned  mass  and  moment  of  inertia. 

As  a  test  was  computed  the  drag  forces  between  the  solid  particles.  The  simulations 
comprise  a  periodic  unit  cell,  2L  x  L  x  L,  with  spheres  located  at  (L  ±  1/2  L,  1/2  L,  1/2  L).  The  two 
spheres  move  with  opposite  velocities  u  and  -u  (thus  there  is  no  momentum  flux). 

For  all  velocities  the  drag  force  was  found  to  be  linear  and  superposable  function  of  u.  In 
Figure  23  the  simulation  results  are  compared  with  integral-equation  solutions  for  an  identical 
geometry,  including  an  exact  calculation  of  the  lubrication  forces.  The  overall  agreement  is  good, 
the  simulation  being  accurate  for  particle  separation  of  1  lattice  space.  Further  on,  the 
hydrodynamic  transport  coefficients  of  equilibrium  distributions  were  calculated  (Ladd  1994a,  b). 
Have  been  computed  the  permeability  of  fixed  random  arrays  of  spheres  (K),  the  collective 
mobility  or  sedimentation  velocity  (pi),  the  short-time  self-diffusion  coefficient  (D$)  and  the 
high-frequency  viscosity  (rj„).  In  Figure  24  the  results  from  the  Lattice-Boltzmann  method 
simulations  are  compared  with  independent  calculations,  based  on  multipole  moment  expansions 
of  the  Oseen  equation.  The  results  when  corrected  for  finite-size  effects  are  in  excellent 
agreement  with  experiment. 

The  self-diffusion  coefficient  Ds  was  also  computed  from  the  average  mean-square 
displacement  or  from  the  velocity  auto  correlation  function.  Figure  25  gives  the  ratio  Ds/D0  for 
various  volume  fractions  and  particle  sizes.  The  agreement  is  quite  good  but  compared  with  the 
steady  non-equilibrium  flows  (figure  24)  larger  particles  are  required  at  higher  densities  for  the 
same  accuracy.  These  discrepancies  could  came  from  the  effects  of  insufficient  shared  boundary 
nodes,  because  for  larger  spheres  the  agreement  is  good.  Also  the  results  for  the  collective 
mobility  are  reported  in  Figure  25  and  the  agreement  is  quite  good. 

These  results  as  well  as  other  reported  recently  (Behrend  1995)  show  that  the  lattice 
Boltzmann  method  as  well  as  lattice  Boltzmann-BGK  is  a  viable  technique  for  quantitative 
simulations  of  hydrodynamicaly  interacting  particles.  Even  using  small  solid  particles,  with  radii 
less  than  5  lattice  spacing,  accurate  results  for  hydrodynamic  transport  coefficients  (permeability, 
sedimentation  velocity,  self-diffusion  coefficient  and  viscosity)  have  been  obtained  over  a  whole 
range  of  packing  fractions. 
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FIGURE  23.  Hydrodynamic  interactions  between  pairs  of  spheres.  The  parallel  and 
perpendicular  friction  coefficients  are  plotted  as  a  function  of  5  =  R/a-2.  The  results  at  5  =  0  are 
for  objects  in  closest  possible  proximity.  The  systems  are  periodic,  with  a  two-sphere  unit  cell. 
The  solid  lines  are  again  solutions  of  the  Stokes  equations  in  the  same  periodic  geometry. 
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FIGURE  24.  Hydrodynamic  transport  coefficients  of  random  arrays  of  spheres.  Results  from 
simulations  of  16  spheres  (with  periodic  boundary  conditions)  are  compared  with  accurate 
numerical  solutions  of  the  Stokes  equations  (Ladd,  1990).  The  lattice-Boltzmann  results  (present 
work)  are  plotted  as  symbols;  results  from  Ladd  (1990)  are  shown  as  solid  lines.  The  statistical 
errors  in  both  sets  of  calculations  are  smaller  than  the  plotting  symbols. 
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FIGURE  25.  Diffusion  coefficients  of  random  arrays  of  spheres.  Self-diffusion  coefficient  Ds 
and  collective  mobility  p  from  simulations  of  16  spheres  (filled  symbols)  and  128  spheres  (open 
symbols)  are  compared  with  numerical  solutions  of  the  Stokes  equations  (Ladd,  1990).  Results 
from  moment  expansions  of  the  Oseen  equation  are  shown  as  solid  lines  (N  =  16)  and  dashed 
lines  (IV  =  128).  The  statistical  errors  in  the  fluctuating  lattice-Boltzmann  simulations  are 
comparable  to  the  size  of  the  plotting  symbols. 


4. 3.  Free  boundary  flows 

The  Lattice-Boltzmann  methods  was  extended  by  a  reservoir  of  rest  particles  to  model 
the  flow  and  pressure  distribution  in  a  lid-driven  cavity  at  high  Reynolds  numbers. 

By  introducing  a  reservoir  of  rest  particles  with  a  density  do,  the  dependence  of  pressure  p 
on  velocity  u  can  be  removed  (Chen  et  al  1992),  The  total  density  of  the  system  is  then  given  by: 

P  =  d0  +  Mcd  (139) 

with  d  of  moving  particles  per  node  and  with  the  usual  LB-BGK  equations  changing  accordingly 
(Miller  1995,  Eggels  1995). 

The  lattice  size  was  150  x  150  mesh  points  was  used  to  study  the  flow  in  a  cavity  with 
different  types  of  boundary  conditions  (uniform  and  non  uniform  Dirichlet,  uniform  and  non 
uniform  Neumann). 

For  Reynolds  number  around  1000  and  a  uniform  shear  flow  applied  at  the  top  of  the 
lattice,  two  strong  vortices  are  formed  on  the  right-hand  side  of  the  cavity  (Figure  26).  For  the 
case  with  the  non  uniform  shear  flow,  only  one  strong  vortices  are  formed  at  the  right  hand  side, 
due  to  a  high  pressure  zone  formation  at  the  end  of  the  cavity. 
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FIGURE  26.  Streamlines  for  Re  =  1000,  aspect  ratio  A  =  10,  and  the  uniform  Dirichlet  boundary 
condition  at  the  top. 
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FIGURE  27.  Streamlines  for  Re  =  1000,  aspect  ratio  A  =  10,  and  the  nonuniform  Dirichlet 
boundary  condition  at  the  top. 
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5.  Conclusions 


The  lattice-gas  models  (assuming  identical  particles  that  hop  from  site  to  site  on  a  regular 
triangular  lattice  and  obeying  simple  collision  rules  that  conserve  mass  and  momentum) 
demonstrated  that  the  full  details  of  real  molecular  dynamics  are  not  necessary  to  create  a 
microscopic  model  with  a  macroscopic  hydrodynamic  behavior. 

Further  on,  the  LGA  became  an  alternative  means  for  the  numerical  simulation  of 
hydrodynamic  flows  and  gave  rise  to  new  ideas  for  constructing  models  of  complex  fluids:  fluid 
mixture  including  interfaces,  exhibiting  phase  transitions  or  multi-phase  flows  through  porous 
media,  etc. 

The  main  advantages  of  LGA’s  techniques  are:  intrinsic  stability,  easy  introduction  of 
boundary  conditions  and  therefore  a  great  facility  to  deal  with  complex  geometries,  simple 
numerical  codes  and  an  intrinsic  parallel  structure  allowing  an  efficient  parallel  implementation. 

The  early  drawbacks  of  the  LGA  models:  statistical  noise  (and  therefore  the  need  for 
spatial/time  averaging),  not-Galilean  invariant,  (the  presence  of  the  g(p)  factor  in  the  nonlinear 
advection  term  of  the  Navier-Stokes  equations),  the  velocity-dependent  pressure  and  the  presence 
of  the  spurious  invariants,  have  been  cured  by  the  recent  developments  of  the  Lattice-Boltzmann 
and  the  Lattice-Boltzmann  BGK  method. 

The  LGA  techniques,  while  developed  initially  for  the  fluid  hydrodynamics,  proved  to  be 
applicable  to  a  wide  range  of  phenomena  through  a  proper  choice  of  the  equilibrium  distribution 
function  (the  actual  choice  being  just  an  example  which  is  suitable  for  the  recovering  of  the 
Navier-Stokes  equations):  reaction-diffusion  systems,  crystal  growth,  heterogeneous  chemical 
reactions,  diffusion,  wave  propagation,  and  magneto-hydrodynamics. 


88 


SELECTED  REFERENCES 


Lattice-gas 

F.  J.  Alexander,  I.  Edrei,  P.L.  Garrido,  J.  L.  Lebowitz  (1992)  Phase  transitions  in  a 
probabilistic  cellular  automaton:  growth  kinetics  and  critical  properties.  J.  Stat.  Phys.  vol.  68,  nos. 
3/4,  pp.  497-514. 

C.  Appert,  D.  Rothman,  S.  Zaleski  (1991)  A  liquid-gas  model  an  a  lattice.  Physica  D,  vol. 
47,  pp. 85-96. 

C.  Appert  (1993)  Gas-liquid  dynamical  phase-change  and  interface  generation  on  lattice 
gases.  Phd.  Thesis,  Ecole  Normale  Sup.,  Univ.  Paris  6,  Paris,  (in  french). 

G.  K.  Batchelor  (1967)  An  introduction  to  fluid  dynamics.  Cambridge  Univ.  Press,  NY. 

D.  Bemardin,  O.  E.  Sero-Guillaume,  C.  H.  Sun  (1991)  Multispecies  2D  lattice  gas  with 
energy  levels:  diffusive  properties.  Physica  D  47,  pp.  169-188. 

B.  M.  Boghosian,  D.  Levermore  (1987)  A  cellular  automaton  for  Burger’  equation.  Complx. 
Syst.  no.  1,  pp.  17-30. 

J.  P.  Boon  (1991)  Statistical  mechanics  and  hydrodynamics  of  lattice  gas  automata:  An 
overview.  Physica  D,  vol.  47,  pp.3-8. 

L.  Brieger,  E.  Bonomi  (1991)  A  stochastic  cellular  automaton  simulation  of  the  non-linear 
diffusion  equation.  Physica  D,  vol.  47,  pp.159-168. 

C.  Burges,  S.  Zaleski  (1987)  Buoyant  mixtures  of  cellular  automaton  gases.  Complx.  Syst. 
vpl.  1,  pp.  31-50. 

C.  K.  Chan,  N.  Y.  Liang  (1990)  Critical  phenomena  in  an  immiscible  lattice-gas  cellular 
automaton.  Europhysics  Letters,  vol.  13,  no.  6,  pp.  495-500. 

H.  Chen,  S.  Chen,  G.  Doolen,  Y,  C.  Lee  (1988)  Simple  lattice  gas  models  for  waves. 
Complx.  Syst.  vol.2,  pp.  259-267. 

S.  Chen,  D.O.  Martinez,  W.  H.  Mattheus,  H.  Chen  (1992)  Magnetohydrodynamics 
computations  with  lattice  gas  automata.  J.  Stat.  Phys.,  vol.  68,  nos.  3/4. 

S.  Chen,  H.  Chen,  G.G.  Doolen,  Y.  C.  Lee,  H.  Rose  (1991)  Lattice  gas  models  for  nonideal 
gas  fluids.  Physica  D,  vol.  47,  pp.  97-1 1 1. 

S.  Chen,  K.  Diemer,  G.  D.  Doolen,  K.  Eggert,  C.  Fu,  S.  Gutman,  K.  Eggert,  B.  J.  Travis 
(1991)  Lattice  gas  automata  for  flow  through  porous  media.  Physica  D,  vol.  47,  pp.  72-84. 

B.  Chopard,  M.  Droz  (1990)  Cellular  automata  model  for  thermo-hydrodynamics,  in 
‘Numerical  methods  for  simulation  of  multi-phase  and  complex  flow”,  Proc.  of  a  Workshop  held  at 
Koninklijeke/Shell-Lab.,  30  May-1  June. 

P.  Clavin,  D.  d’Humieres.  P.  Lallemand,  Y.  Pomeau  (1990)  Cellular  automata  for 
hydrodynamics  with  free  boundarries  in  two  and  three  dimensions,  in  “Lattice  gas  methods  for 
partial  differential  equations”,  Eds.  Dolen  et  al.,  Addison-Weslesy  Publ.  Co. 

89 


K.  A.  Cliffe,  R.D.  Kingdon,  P.  Schofield,  P.J.  Stopford  (1991)  Lattice  gas  simulation  of  free 
boundary  flows.  Physica  D,  vol.  47,  pp.  275-280. 

A.  Clouqueur,  D.  d’Humieres  (1987)  RAP1,  a  cellular  automaton  machine  for  fluid 
dynamics.  Complx.  Syst.  vol.l,  pp.  585-597. 

R.  Comubert,  D.  d’Humieres,  D.  Levermore  (1991)  A  Knudsen  layer  theory  for  lattice  gases. 
Physica  D,  vol.  47,  pp.  241-259. 

B.  Dubrulle,  U.  Frisch,  M.  Henon,  J-P.  Rivet  (1991)  Low- viscosity  latice  gases.  Physica  D, 
vol.  47,  pp.  27-29. 

M.  H.  Ernst,  S.  P.  Das,  D.  Kammann  (1990)  Heat  conductivity  in  a  thermal  LGCA.  in 
‘Numerical  methods  for  simulation  of  multi-phase  and  complex  flow”,  Proc.  of  a  Workshop  held  at 
Koninklijeke/Shell-Lab.,  30  May-1  June. 

U.  Frisch,  J.  P.  Rivet  (1986)  Lattice  gas  hydrodinamics,  Green-Kubo  formula.  Comptes 
Rendus  303,  II,  pp.  1065. 

U.  Frisch,  B.  Hasslacher,  Y.  Pomeau  (1986)  Lattice-gas  automata  for  the  Navier-Stokes 
equation.  Phys.  Rev.  Lett.,  vol.  56,  no.  14,  pp.  1505-1508. 

U.  Frisch,  D.  d’Humieres,  B.  Hasslacher,  P.  Lallemand,  Y.  Pomeau,  J.-P.  Rivet  (1987) 
Lattice  gas  hydrodynamics  in  two  and  three  dimensions.  Complx.  Syst.  no.l,  pp.  649-707. 

R.  Gatignol  (1987)  The  hydrodynamical  description  for  a  dicrete  velocity  model  of  gas. 
Complx.  Syst.  no.l,  pp.  709-725. 

A.  K.  Gunstensen,  D.  H.  Rothman  (1991)  A  Galilean-invariant  immiscible  lattice  gas. 
Physica  D,  vol.  47,  pp.  53-63. 

R.  Gutfraind,  A.  Hansen  (1994)  Study  of  hydrodynamic  permeability  of  fractures  using 
lattice  gas  automata.  Techn.  Rap.  Univ.  of  Trondheim,  UNIT,  no.3. 

R.  Gutfraind,  I.  Ippolito  (1995)  Study  of  tracer  dispersion  in  self-affine  fractures  using 
lattice-gas  automata.  Phys.  Fluids,  vol.7,  no. 8,  pp.  1938-1948. 

B.  Hasslacher  (1987)  Discrete  fluids.  Los  Alamos  Science  Special  Issue  (preprint). 

T.  Hatori,  D.  Montgomery  (1987)  Transport  coefficients  for  magnetohydrodinamic  cellular 
automata.  Complex  Syst.  vol.  1,  pp.735-752. 

F.  Hayot  (1987)  The  effect  of  Galilean  non-invariance  in  lattice  gas  automaton  one¬ 
dimensional  flow.  Complx.  Syst.,  no.  1,  pp.  753-761. 

F.  Hayot  (1991)  Fingering  instability  in  a  lattice  gas.  Physica  D,  vol  47,  pp.  64-71. 

M.  Henon  (1987a)  Viscosity  of  a  lattice  gas.  Complx.  Syst.  no.l,  pp.  763-789. 

M.  Henon  (1987b)  Isometric  collision  rules  for  hte  four-dimensional  FCHC  lattice  gas. 
Compl.  Syst.  vol.  1,  pp.  475-494. 

M.  Henon  (1992)  Implementation  of  the  FCHC  lattice  gas  model  on  the  connection  machine. 
J.  Stat.  Ohys.,  vol.  68,  nos.  3/4,  pp.  353-377. 


90 


F.  Higuera  (1990)  Modified  lattice  gas  method  for  high  Reynolds  number  incompressible 
flows,  in  “Numerical  methods  for  simulation  of  multi-phase  and  complex  flow”,  Proc.  of  a 
Workshop  held  at  Koninklijeke/Shell-Lab.,  30  May-1  June. 

M.  A.  van  der  Hoef,  D.  Frenkel  (1991)  Tagged  particle  diffusion  in  3D  lattice  gas  automata. 
Physica  D,  vol.  47,  pp.  191-197. 

D.  d’Humieres,  P.  Lallemand  (1987)  Numerical  simulations  of  hydrodynamics  with  lattice 
gas  automata  in  two  dimensions.  Complx.  Syst.  no.  1,  pp.  599-632. 

D.  d’Humieres,  P.  Lallemand,  G.  Searby  (1987)  Numerical  experiments  on  lattice  gases: 
Mixtures  and  galilean  invariance.  Complx.  Syst.  no.  1,  pp.  633-647. 

D.  d’Humieres,  P.  Lallemand  (1990)  Flow  of  a  lattice  gas  between  two  parallel  plates: 
Development  of  the  Poiseuille  profile,  in  “Lattice  gas  methods  for  partial  differential  equations”, 
Eds.  Dolen  et  al.,  Addison- Weslesy  Publ.  Co. 

D.  Jeulin  (1990)  Flow  and  diffusion  in  random  porous  media  from  lattice  gas  simulations.in 
“Numerical  methods  for  simulation  of  multi-phase  and  complex  flow”,  Proc.  of  a  Workshop  held  at 
Koninklijeke/Shell-Lab.,  30  May-1  June. 

L.  P.  Kadanoff,  G.  R.  McNamara,  G.  Zanetti  (1987)  A  Poiseuille  viscometer  for  lattice  gas 
automata.  Complx.  Syst.  no.l,  pp.  791-803. 

J.  M.  V.  A.  Koelman,  M.  Nepveu  (1990)  Darcy  flow  in  porous  media:  Cellular  automata 
simulations,  in  "Numerical  methods  for  simulation  of  multi-phase  and  complex  flow",  Proc.  of  a 
Workshop  held  at  Koninklijeke/Shell-Lab.,  30  May-1  June. 

G.  A.  Kohring  (1992)  The  cellular  automata  approach  to  simulating  fluid  flows  in  porous 
media.  Physica  A,  vol.  186,  pp.  97-108. 

A.  W.  Kolmogorov  (1941)  C.  R.  Acad.  Sci.  USSR,  30,  pp.  301-538. 

R.  H.  Kraichnan  (1967)  Phys.  Fluids,  10,  pp.  1417. 

L.  D.  Landau  and  E.  M.  Lifschitz  (1959)  Fluid  mechanics.  Pergamon  Press.  NY. 

F.  Liu,  N.  Goldenfeld  (1991)  Deterministic  lattice  models  for  diffusion-controlled  crystal 
growth.  Physica  D,  vol.  47,  pp.  124-131. 

J.  F.  Lutsko,  J.  P.  Boon,  J.  A.  Somers  (1990)  Lattice  gas  automata  simulations  of  viscous 
fingering  in  porous  media,  in  "Numerical  methods  for  simulation  of  multi-phase  and  complex 
flow",  Proc.  of  a  Workshop  held  at  Koninklijeke/Shell-Lab.,  30  May-1  June. 

P.  Manneville,  N.  Boccara,  G.Y.  Vichniac,  R.  Bidaux,  eds.  (1989)  Cellular  automata  and 
modeling  of  complex  physical  systems.  Springer-Verlag. 

N.  Margolus,  T.  Toffoli  (1990)  Cellular  automata  machines,  in  “Lattice  gas  methods  for 
partial  differential  equations”,  eds.  Doolen  et  al.,  Addison-wesley  Publ.  Co. 

G.  R.  McNamara  (1990)  Diffusion  in  a  lattice  automaton.  Europhysics  Letters,  vol.  12,  no.4, 
pp.  329-334. 

P.  Mora  (1992)  The  lattice  Boltzmann  phononic  lattice  solid.  J.  stat.  Phys.,  vol.  68,  nos.  3/4. 


91 


F.  R.  Mujica-Femandez  (1991)  Lattice  gas  wind  tunnel.  M.S.  Thesis,  Dep.  Nucl.  Eng.,  MIT, 
Cambridge,  MA. 

T.  Naitoh,  M.H.  Ernst,  M.A.  van  der  Hoef  (1990)  Effects  of  sound  modes  on  the  VACF  in 
cellular  automaton  fluids,  in  ‘Numerical  methods  for  simulation  of  multi-phase  and  complex  flow”, 
Proc.  of  a  Workshop  held  at  Koninklijeke/Shell-Lab.,  30  May-1  June. 

J.  von  Neumann  (1966)  Theory  of  self-reproducting  automata.  (Univ.  of  Ilinois  press), 

J.  F.  Olson  (1995)  Two-fluid  flow  in  sedimentary  rock:  complexity,  transport  and  simulation. 
PhD  Thesis,  Dep.  of  Earth,  Atms.  and  Planetary  Sciences,  MIT,  Cambridge,  MA. 

S.  A.  Orszag,  V.  Yakhot  (1986)  Reynolds  number  scaling- of  cellular  automaton 
hydrodynamics.  Phys.  Rev.  Lett.,  vol.  56,  no.  16,  pp.  1691-1693. 

V.  Pot  (1994)  Microscopical  study  of  flow  and  phase-change  in  porous  media  via  lattice  gas 
models.  Phd  Thesis,  Univ.  Paris  6,  (in  french). 

Y.  H.  Qian,  D.  d’Humieres,  P.  Lallemand  (1992)  Diffusion  simulation  with  a  deterministic 
one-dimensional  lattice-gas  model.  J.  Stat.  Phys.,  vol.  68,  nos.  3/4,  pp.  563-573. 

F.  Reif  (1988)  Fundamental  of  statistical  and  thermal  physics.  Mc.Graw-Hill,  NY. 

J.  -P.  Rivet,  U.  Frisch  (1987)  Green-Kubo  formalism  for  lattice  gas  hydrohynamics  and 
Monte-Carlo  evaluation  of  shear  viscosities.  Complx.  Syst.  no.  1,  pp.  839-851. 

J.-P.  Rivet,  M.  Henon,  U.  Frisch,  D.  d’Humieres  (1988)  Simulating  fully  three-dimensional 
external  flow  by  lattice  gas  methods.  Europhysiscs  Letters,  vol.  7,  no.  3,  pp.  231-236. 

D.  Rothman,  J.  M.  Keller  (1988)  Immiscible  cellular-automaton  fluids.  J.  Stat.  Phys.,  vol.  52, 
nos.  3/4,  pp.  1119-1127. 

D.  Rothman  (1988)  Cellular-automaton  fluids:  A  model  for  flow  in  porous  media. 
Geophysics,  vol.  53,  no.4,  pp.  509-518. 

D.  Rothman  (1990)  Macroscopic  laws  for  immiscible  two-phase  flow  in  porous  media: 
Results  from  numerical  experiments.  J.  Geoph.  Res.,  vol.  95,  no.  B6,  pp.  8663-8674. 

D.  Rothman  (1992)  Simple  models  of  complex  fluids,  in  “Microscopic  simulations  of 
complex  hydrodynamics”,  eds.  M.  Marechal  and  B.  Holian,  Plenum  Press. 

D.  Rothman,  L.  P.  Kadanoff  (1994)  Bubble,  bubble,  boil,  and  trouble.  Comp.  Phys.,  vol.  8, 
no.  2,  pp.  199-204. 

D.  Rothman,  S.  Zaleski  (1994)  Lattice  gas  models  of  phase  separation:  interfaces,  phase 
transitions,  and  multiphase  flow.Rev.  Mod.  Phys.,  vol.  66,  pp.  1417-1479. 

D.  Rothman  and  S.  Zaleski  (1966)  Lattice-gas  cellular  automata:  Simple  models  of  complex 
hydrodynamics,  (in  print). 

T.  Shimomura,  G.  D.  Doolen,  B.  Hasslacher,  C.  Fu  (1987)  Calculation  using  lattice  gas 
techniques.  Los  Alamos  Sci.,  Special  Issue,  pp.  210-210. 

J.  A.  Somers,  P.  C.  Rem  (1990)  Obtaining  numerical  results  from  the  3D-FCHC  lattice 
gas,  in  "Numerical  methods  for  simulation  of  multi-phase  and  complex  flow",  Proc.  of  a 
Workshop  held  at  Koninklijeke/Shell-Lab.,  30  May-1  June. 

92 


J.  A.  Sommers,  P.  C.  Rem  (1991)  Flow  computation  with  lattice  gases,  in  R,  V.  Oliemans 
(ed.)  Comp.  Fluid  Dyn.,  Kluwer  Ac.  Publ. 

J.  A.  Sommers,  P.  C.  Rem  (1991)  Analysis  of  surface  tension  in  two-phase  lattice  gases. 
Physica  D,  vol.  47,  pp.  39-46. 

G.  Y.  Vichniac  (1990)  cellular  automata  and  complex  systems,  in  “Microscopic  simulations 
of  complex  hydrodynamics”,  eds.  M.  Marechal  and  B.  Holian,  Plenum  Press. 

G.  Zanetti  (1990)  Lattice  gas  automata:  Comparison  of  simulation  and  theory,  in 
“Microscopic  simulations  of  complex  hydrodynamics”,  eds.  M.  Marechal  and  B.  Holian,  Plenum 
Press. 

G.  Zanetti  (1991)  Counting  hydrodynamic  modes  in  lattice  gas  automata.  Physica  D,  vol.  47, 
pp.  30-35. 

J.T.  Wells,  D.R.  Janecky,  BJ.  Travis  (1991)  A  lattice  gas  automata  model  for  heterogenous 
chemical  reactions  at  mineral  surfaces  and  in  pore  networks,  Physica  D,  fol.  47,  pp.l  15-123. 

S.  Wolfram  (1986)  Cellular  automaton  fluids  1:  Basic  Theory.  J.  Stat.  Phys.,  vol  45,  nos.3/4, 
pp.  471-526. 

Lattice-gas  Boltzmann 

R.  Benzi,  S.  Succi,  M.  Vergassola  (1992)  The  lattice  Boltzmann  equation:  theory  and 
applications.  Phys.  Rep.  222,  no.  3,  pp.  145-197,  North-Holland. 

J.G.M.  Eggels,  J.A.  Sommers  (1995)  Numerical  simulation  of  free  convective  flow  using  the 
lattice-Boltzmann  scheme.  Int.  J.  Heat  and  Fluid  Flow  (pre-print). 

J.  G.  M.  Eggels  (1995)  Direct  and  large-eddy  simulation  of  turbulent  fluid  flow  using  the 
lattice-Boltzmann  scheme.  Int.  J.  Heat  and  Fluid  Flow  (pre-print). 

U.  Frisch  (1991)  Relation  between  the  lattice  Boltzmann  equation  and  the  Navier-Stokes 
equations.  Physica  D,  vol.  47,  pp.  231-232. 

B.  Ferreol,  D.  Rothman  (1995)  Lattice-Boltzmann  simulations  of  flow  through  Fontainebleu 
sandstone.  Trans,  in  Porous  Media,  (pre-print). 

I.  Ginzbourg.P.M.  Adler  (1994)  Boundary  flow  condition  analysis  for  the  three-dimensional 
lattice-Boltzmann  model.  J.  Phys.  II,  France,  vol.  4,  pp.  191-214. 

A.  K.  Gunstensen  (1992)  Lattice-Boltzmann  studies  of  multi-phase  flow  through  porous 
porous  media.  Phd  Thesis,  Dep.  of  Earth  Atms.  and  Planetary  Science,  MIT,  Cambridge,  MA. 

A.  K.  Gunstensen,  D.H.  Rothman  (1993)  Lattice-Boltzmann  studies  of  immiscible  two-phase 
flow  through  porous  media.  J.  Geoph.  Res.,  vol.  98,  no.  B4,  pp.  6431-6441. 

F.J.  Higuera,  S.  Succi,  R.  Benzi  (1989)  Lattice  gas  with  enhanced  collisins.  Europhys.  Lett. 
vol.9,  no.4,pp.  345-349. 

F.  J.  Higuera,  J.  Jimenez  (1989)  Boltzmann  approach  to  lattice  gas  simulations.  Europhys. 
Lett.,  vol.9,  no.7,  pp.  663-668. 


93 


F. J.  Higuera,  S.  Sued  (1989)  Simulating  the  flow  around  a  circular  cylinder  with  a  lattice 
Boltzmann  equation.  Europhysics  Letters,  vol.  8,  no.  6,  pp.  517-521. 

R.  Holme,  D.H.  Rothman  (1992)  Lattice-gas  and  Lattice-Boltzmann  models  of  miscible 
fluids.  J.  of  Stat.  Physics,  vol.  68,  nos.  3/4,  pp.  409-429. 

D.  d’Humieres  (1992)  Generalized  lattice  Boltzmann  equations,  (pre  print) 

G.  R.  McNamara,  G.  Zanetti  (1988)  Use  of  the  Boltzmann  equation  to  simulate  lattice-gas 
automata.  Phys.  Rev.  Lett.,  vol.61,no.  20. 

D.R.  Noble,  S.  Chen,  J.  G.  Georgiadis,  R.  O.  Buckius  (1995)  A  consistent  hydrodynamic 
boundary  condition  for  the  lattice  Boltzmann  method.  Phys.  of  Fluids,  vol.  -7,  no.l,  pp.  203-207. 

U.  d’Ortona,  D.  Salin,  M.  Cieplak,  J.R.  Banavar  (1995)  Interfacial  phenomena  in  Boltzmann 
cellular  automata.  Europhys.  Lett.,  vol.  28,  no.5,  pp.  317-322. 

U.  d’Ortona,  D.  Salin,  M.  Cieplak,  R.  Rybka,  J.R.  Banavar  (1995)Two-color  nonlinear 
Boltzmann  cellular  automata:  Surface  tension  and  wetting.  Phys.  rev.  E,  vol.  51,  no.4,  pp.  3718- 
3728. 

R. B.  Rybka,  M.  Cieplak,  D.  Salin  (1995)  Boltzmann  cellular  automata  studies  of  the 
spinodal  decomposition.  Physica  A,  vol.  222,  pp.  105-1 18. 

J.P.  Rivet,  U.  Frisch  (1986)  Lattice  gas  automata  in  the  Boltzmann  approximation.  Compte 
Rendus  302,  n,  pp.  267. 

P.  A.  Skordos  (1995)  Modeling  flue  pipes:  Subsonic  flow,  lattice  Boltzmann  and  parallel 
distributed  computers.  Phd  Thesis,  Dep.  Elec.  Eng.  and  Comp.  Sci.,  MIT,  Cambridge,  MA. 

S.  Succi,  R.  Benzi,  F.  Higuera  (1991)  The  lattice  Boltzmann  equation:  A  new  tool  for 
computational  fluid-dynamics.  Physca  D,  vol.  47,  pp.219-230. 

C.M.  Teixeira  (1992)  Continuum  limit  of  lattice  fluid  dynamics.  Phd  Thesis,  Dep.  Nucl. 
Eng.,  MIT,  Cambridge,  MA. 


Lattice-gas  BGK 

F.  J.  Alexander,  H.  Chen,  S.  Chen,  G.D.  Doolen  (1992)  Lattice  Boltzmann  model  for 
compressible  fluids.  Physical  Rev.  A,  vol.  46,  no.  4,  pp.  1967-1970. 

O.  Behrend,  R.  Harris,  P.  B.  Warren  (1994)  Hydrodynamic  behavior  of  lattice  Bhatnagar- 
Gross-Krook  models.  Phys.  Rev.  E,  vol.  50,  no.6,  pp.  4586-4595. 

O.  Behrend  (1995)  Solid-fluid  boundaries  in  particle  suspension  simulations  via  the  lattice 
Boltzmann  method.  Phys.  Rev.  E,  vol. 52,  no.  1,  pp.l  164-1 175. 

B.  M.  Boghosian,  W.  Taylor  (1995)  Correlations  and  renormalizations  in  lattice  gases.  Phys. 
Rev.  E,  vol.  52,  no.l,  pp.  510-554. 

H.  Chen,  S.  Chen,  W.H.  Mattheus  (1992)  Recovery  of  the  Navier-Stokes  equations  using  a 
lattice-Boltzmann  method.  Phys.  Rev.  A,  vol.  45,  no.  8,  pp.  5339-5342. 

94 


S.  Chen,  Z.  Wang,  X.  Shan,  G.D.  Doolen  (1992)  Lattice  Boltzmann  computational  fluid 
dynamics  in  three  dimensions.  J.  of  Stat.  Phys.,  vol.  68,  nos.  3/4. 

Y.  Chen,  H.  Ohashi,  M.  Akiyama  (1994)  Thermal  lattice  Bhatnagar-Gross-Krook  model 
without  nonlinear  deviations  in  macrodynamic  equations.  Phys.  Rev.  E,  vol.  50,  no.4,  pp.  2776-2783 

E.G.  Flekkoy,  J.  Feder,  T.  Jossang  (1992)  Lattice  gas  simulations  of  osmosis.  J.  of  Stat. 
Phys.,  vol.  68,  nos.  3/4. 

E.G.  Flekkoy  (1993)  Lattice  BGK  for  miscible  fluids,  (preprint). 

E.G.  Flekkoy  (1993)  Modelling  miscible  fluids.  Phd  Thesis,  Centre  for  Advance  Study, 
Oslo,  Norway. 

M.  Gerits,  M.H.  Ernst  (1993)  Lattice-gas  automata  with  attractive  and  repulsive  interactions. 
Phys.  Rev.  E,  vol.  48,  no.2,  pp.  988-999. 

S.  Hou  (1995)  Lattice  Boltzmann  method  for  incompressible  viscous  flow.  Ph.D..  Thesis, 
Kansas  State  Univ.,  Manhattan,  Kansas. 

J.  M.  V.  A.  Koelman  (1991)  A  simple  lattice  Boltzmann  scheme  for  Navier-Stokes  fluid 
flow.  Europhys.  Lett.,  vol.  15,  no.6,  pp.  603-607. 

D.  O.  Martinez,  S.  Chen,  W.H.  Matthaeus  (1993)Lattice  Boltzmann  magneto 
hydrodynamics.  Phys.  Plasmas,  vol.l,  no.  6,  pp.  1850-1867. 

W.  Miller  (1995)  Flow  in  the  driven  cavity  calculated  by  the  lattice  Boltzmann  method. 
Phys.  Rev.  E,  vol.51,  no.4,  pp.  3659-3669. 

A.C.  Ladd  (1990)  Hydrodynamic  interactions  and  transport  coefficients  in  a  suspension  of 
spherical  particles,  Microsc..  Sim.  of  Complex.  Flows,  (ed.  M.  Mareschal),  NY. 

A.C.  Ladd  (1991)  Dissipative  and  fluctuating  hydrodynamic  interactions  between  suspended 
solid  particles  via  Lattice-gas  Cellular  Automata,  in  M.Meyer  and  V.  Pontikis  ed.,  Computer  Sim.  in 
mat.  Sci.,  Kluwer  Acad.  Pub.,  pp.  481-504. 

A.  J.  C.  Ladd  (1993)  Short-time  motion  of  colloidal  particles:  Numerical  simulation  via  a 
fluctuating  lattice-Boltzmann  equation.  Phys.  Rev.  Lett.,  vol.  70,  no.  9,  pp.  1339-1342. 

A.C.  Ladd  (1994,  a)  Numerical  simulations  of  particulate  suspensions  via  a  discretized 
Boltzmann  equation.  Part  1.  Theoretical  foundation.  J.  Fluid  Mech.,  vol.  221,  pp.  285-309, 
Cambridge  Univ.  Press. 

A.C.  Ladd  (1994,  b)  Numerical  simulations  of  particulate  suspensions  via  a  discretized 
Boltzmann  equation.  Part  2.  Numerical  Results.  J.  Fluid  Mech.,  vol.  221,  pp.  285-309,  Cambridge 
Univ.  Press. 

A.  C.  Ladd,  Hu  Gang,  J.X.  Zhu,  D.  A.  Weitz  (1995)  Time-dependent  collective  diffusion  of 
colloidal  particles.  Phys.  Rev.  Lett.,  vol.74,  no.  2,  pp.  318-321. 

A.  Lawniczack,  D.  Dab,  R.  Kapral,  J.P.  Boon  (1991)  Reactive  lattice  gas  automata.  Physica 
D,  vol.47,  pp.  132-158. 

D.  R.  Noble,  S.  Chen,  J.  G.  Georgiadis,  R.  O.  Buckius  (1995)  A  consistent  hydrodynamic 
boundary  condition  for  the  lattice  Boltzmann  method.  Phys.  Fluids,  vol.  7,  no.  1,  pp.  203-209. 

95 


Y.  H.  Qian,  D.  d’Humieres,  P.  Lallemand  (1992)  Lattice  BGK  for  Navier-Stokes  equations. 
Europhys.  Lett,  vol.17,  no.  6,  pp.  479-484. 

Y.  H.  Qian,  S.  A.  Orszag  (1993)  Lattice  BGK  models  for  the  Navier-Stokes  equation: 
Nonlinear  deviation  in  compressible  regimes.  Europhys.  Lett.,  vol.  21,  no.  3,  pp.  255-259. 

S.  Schwarzer  (1995)  Sedimentation  and  flow  through  porous  media:  Simulating  dynamically 
coupled  discrete  and  continuum  phases.  Phys.  Rev.  E,  vol.  52,  no.6,  pp.  6461-6475. 

X.  Shan,  H.Chen  (1992)  Lattice  Boltzmann  model  for  simulating  flows  with  multiple  phases 
and  components.  Phys..  Rev.  E,  vol.  ,  no.  ,  pp.  1815-1819. 

X.  Shan,  H.  Chen  (1995)  Simulation  of  non  ideal  gases  and  liquid-gas  phase  transitions  by 
the  lattice  Boltzmann  equation.  Phys.  Rev.  E,  no.  49,  vol.4,  pp.  2941-2948. 

P.A.  Skordos  (1993)  Initial  and  boundary  conditions  for  the  lattice  Boltzmann  method.  Phys. 
Rev.  E,  vol.  48,  no.6,  pp.  4823-4842. 


96 


Data  Report  for  the  1995  Wind  River  Mountains  - 
Green  River  Basin,  Wyoming,  Seismic  Refraction 
Profile 


Stephan  H.  Koester 
Boston  College 
Phillips  Laboratory/GPE 
29  Randolph  Road 
Hanscom  AFB,  MA  01731 


John  J.  Cipar 
Phillips  Laboratory/GPE 
29  Randolph  Road 
Hanscom  AFB,  MA  01731 


14  May  1996 


97 


The  1995  Wind  River  Mountains  -  Green  River 
Basin,  Wyoming,  Seismic  Refraction  Profile 


1.0  Introduction 

During  the  summer  of  1995,  the  Earth  Sciences  Division  of  Phillips 
Laboratory  (PL/GPE)  undertook  an  extensive  seismic  experiment  in 
southwestern  Wyoming.  The  experiment  had  two  parts:  (a)  recording  a 
seismic  refraction  profile  across  the  Green  River  Basin  using  explosions 
fired  east  of  Lander,  Wyoming  as  part  of  the  Deep  Probe  experiment;  and 
(b)  installation  and  operation  of  a  large-aperture  array  during  August  and 
September.  A  description  of  the  technical  aspects  of  the  refraction  profile 
is  the  subject  of  this  first  report;  a  separate  report  covers  data  recorded  by 
the  large-aperture  array. 

In  August,  1995,  several  US  and  Canadian  universities  along  with  the 
Canadian  Geological  Survey  performed  the  Deep  Probe  Experiment,  an 
ultra-large  scale  active  seismic  experiment  in  western  North  America 
(Henstock  et  al,  1995).  The  main  Deep  Probe  profile  was  oriented  north- 
south  from  Edmonton,  Alberta,  to  Crownpoint,  New  Mexico,  a  distance  of 
1900  km.  An  intermediate  shot  point  was  located  approximately  50  km 
east  of  Lander,  Wyoming,  and  provided  the  sources  for  the  Wind  River 
Mountains  -  Green  River  Basin  Seismic  Refraction  Profile  described  in  this 
report.  This  profile  consisted  of  47  shot-station  points  extending  from  Big 
Sandy,  Wyoming,  west  to  the  Idaho-Wyoming  border,  a  distance  of  about 
150  km. 

We  plan  to  use  the  data  from  this  experiment  to  study: 

•  regional  wave  propagation  in  the  central  Rocky  Mountains 

•  spatial  variability  of  earthquake/explosion  discriminants 

•  crustal  structure  beneath  the  Wind  River  Mountains  and  Green 
River  Basin 

2.0  Previous  Geophysical  Work 

The  Green  River  Basin  is  an  extensive  Cenozoic  sedimentary  basin 
bounded  by  the  Precambrian  Wind  River  Mountains  on  the  east,  the 
Wyoming  Range  on  the  west,  and  the  Uinta  Mountains  of  Utah  to  the  south. 
The  southern  part  of  the  basin  contains  the  Mesozoic  Rock  Springs  uplift. 
While  the  underlying  rocks  are  of  various  ages,  the  overall  structure  of 
mountain  bounded  basin  was  formed  during  the  late  Cretaceous-early 
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Tertiary  Laramide  orogeny.  The  Wind  River  Mountains  are  a  thrust- 
faulted  basement  anticline  that  over  rode  the  eastern  part  of  the  Green 
River  basin  (Smithson  et  al,  1979). 

The  Green  River  basin  abuts  the  eastern  edge  of  the  Intermountain 
Seismic  Belt  with  considerable  seismicity  to  the  south  and  west  of  the  basin 
(Pechmann  et  al,  1995).  Previous  geophysical  studies  include  a  seismic 
refraction  profile  from  American  Falls  Reservoir,  Idaho,  to  Flaming 
Gorge  Reservoir,  Wyoming,  with  an  intermediate  shot  point  at  Bear  Lake, 
Idaho-Utah  (Prohdel,  1979).  Pechmann  et  al  (1995)  use  the  Prohdel 
(1979)  P-wave  velocity-depth  model  to  infer  S-wave  velocities  and 
densities.  This  model  has  crustal  thickness  of  40  km  and  is  underlain  by 
7.9  km/s  mantle  material.  Braile  et  al  (1974)  interpreted  a  single-ended 
refraction  profile  extending  from  the  Bingham  Canyon  copper  mine  near 
Salt  Lake  City,  Utah,  across  the  Green  River  Basin  to  the  Wind  River 
Mountains.  They  infer  that  the  crustal  thickness  is  40  km  or  greater 
beneath  the  Green  River  Basin  and  southern  Rocky  Mountains.  Smithson  et 
al  (1979)  and  Brewer  et  al  (1980)  discuss  COCORP  deep  seismic 
reflection  data  collected  across  the  southern  end  of  the  Wind  River 
Mountains  and  adjacent  Green  River  and  Wind  River  Basins.  These 
observations  indicate  the  shallow  overthrust  nature  of  the  Wind  River 
Mountains. 

3.0  The  Wind  River  Mountains  -  Green  River  Basin  Seismic 

Refraction  Profile 

The  Wind  River  Mountains  -  Green  River  Basin  Seismic  Refraction 
Profile  was  deployed  during  the  Deep  Probe  project  of  August  1995  to 
image  the  crustal  structure  beneath  the  Wind  River  Mountains  and  the 
Green  River  Basin.  Unlike  the  Deep  Probe  profiles,  the  Wind  River 
Mountains  -  Green  River  Basin  Seismic  Refraction  Profile  was  oriented 
east-west.  The  refraction  profile  begins  on  the  west  side  of  the  Wind  River 
Mountains,  135  km  from  the  quarry,  and  extends  in  the  general  azimuth 
of  265°  across  the  Green  River  Basin  and  the  Wyoming  Range  to  a  total 
distance  of  about  280  km  from  the  quarry  (Figure  1).  A  total  of  47 
stations  were  installed,  spaced  approximately  3.2  km  apart  from  each 
other. 


The  refraction  study  utilized  two  Deep  Probe  explosions  east  of 
Lander,  Wyoming,  as  its  energy  source.  The  first  explosion  (labeled  143) 
was  about  15,000  lbs.  of  ammonium  nitrate/fuel  oil  (ANFO)  detonated  on 
the  bottom  of  a  water  filled  quarry  (G.  R.  Keller  personal  communication). 
The  second  explosion  (labeled  243)  was  fired  one  week  later  with 
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approximately  the  same  amount  of  ANFO  in  the  same  location  (Table  1). 
The  large  source  size  and  excellent  coupling  in  the  quarry  provided 
exceptional  signal-to-noise  ratios  across  the  profile. 

3.1  Installation  and  Instrumentation  of  the  Profile 

Three  different  types  of  portable  stations  were  deployed  along  the 
profile  (Table  2).  Figure  2  shows  the  basic  setup  of  a  refraction  profile 
station.  In  particular,  the  station  setup  displayed  uses  RefTek  and 
GeoSpace  instruments,  but  can  also  be  used  as  a  general  guideline  for 
every  station  configuration  used  in  the  profile.  Common  to  each 
configuration  is  a  seismometer  buried  about  one  foot  deep,  a  data 
acquisition  recorder,  a  power  supply,  and  a  timing  system. 

The  first  configuration  consisted  of  3  component,  1.0  Hz  GeoSpace 
HS-10-lb  seismometers  recorded  and  digitized  by  a  24  bit  Refraction 
Technology  data  acquisition  system  (DAS).  The  DAS  contained  a  hard  disk 
drive  to  store  the  digitized  signal  and  was  powered  by  two  12v  gel  cells.  A 
global  positioning  satellite  (GPS)  clock  was  connected  to  the  DAS  to 
provide  accurate  timing  and  geographic  location.  The  DAS  was  set  up  to 
record  a  continuous  single  data  stream  at  100  samples  per  second  consisting 
of  3  data  channels  at  24  bit  resolution  with  a  preamplifier  gain  of  32.  The 
nominal  sensitivity  of  the  HS-10-lb  seismometer  is  600  v/m/sec. 

The  second  station  configuration  deployed  along  the  profile  included 
a  single  vertical  component  1.0  Hz  GeoSpace  HS- 10/lb  seismometer 
digitized  by  a  Terra  Technology  recorder  with  WWVB  radio  timing.  The 
Terra  Technology  recorders  used  an  alarm  clock  to  start  data  acquisition 
before  the  anticipated  chemical  explosion.  Once  triggered,  the  recorders 
acquired  data  for  about  10  to  15  minutes  and  stored  the  data  onto  a  cassette 
tape.  The  Terra  Technology  recorders  were  set  to  record  one  data  channel 
with  12  bit  resolution  using  a  static  gain  of  either  100  or  1000  at  ,100 
samples  per  second. 

The  final  station  configuration  consisted  of  1.0  Hz  MARK-L-4C-3D 
geophones  recorded  and  digitized  by  a  Teledyne  Brown  Engineering 
Portable  Data  Acquisition  System  (PDAS).  The  PDAS  used  GPS  for  both 
timing  and  geographic  location.  The  PDAS  system  was  configured  to 
record  for  15  minutes  starting  at  the  shot  origin  time.  The  parameters 
were  set  to  acquire  3  data  channels  at  14/2  bit  resolution  with  a 
preamplifier  gain  of  1  at  100  samples  per  second.  The  PDAS  data  was 
later  resampled  at  125  samples  per  second  to  conform  with  other  Deep 
Probe  data  sets. 
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In  addition,  data  from  three  large  aperture  array  stations  are 
included  in  the  data  set.  These  stations  are  spaced  approximately  50  km 
from  each  other  in  nearly  a  straight  line  and  are  located  in  the  same 
azimuth  as  the  refraction  profile. 

3.2  Operation  of  Profile 

The  profile  was  deployed  in  two  stages  to  increase  the  offset  from 
the  shot  point.  Since  the  two  chemical  explosions  were  at  the  same 
location,  all  of  the  instruments  were  re-deployed  to  different  locations 
along  the  profile  resulting  in  a  longer  and  denser  profile.  The  first 
deployment  occurred  on  August  8th,  one  day  before  the  first  ANFO 
detonation.  Fourteen  stations  were  installed  for  this  deployment.  Seven 
stations  were  equipped  with  RefTek  recorder  configurations  and  seven 
stations  were  equipped  Terra  Technology  recorder  configurations.  The 
deployment  of  stations  began  at  the  entrance  to  the  Bridger  National  Forest 
with  a  Refraction  Technology  recorder  configuration,  continued  through 
the  national  forest  alternating  recorder  configurations,  and  ended  at  the 
border  of  Wyoming  and  Idaho  with  a  Terra  Technology  recorder 
configuration.  The  stations  were  removed  the  following  day  (August  9th) 
in  the  reverse  order  in  which  they  were  deployed.  Thus,  the  stations  that 
were  retrieved  first  include  less  data  than  those  retrieved  last. 

The  second  deployment  occurred  on  August  16th,  one  day  before  the 
second  explosion.  Thirty  three  stations  were  deployed  including  Refraction 
Technology,  Terra  Technology  and  PDAS  recording  configurations.  The 
second  deployment  stretched  from  the  edge  of  the  Wind  River  Mountains 
west  to  the  entrance  of  the  Bridger  National  Forest.  The  stations  were 
removed  the  following  day  on  August  17th. 

Stations  were  located  along  existing  roads  which  necessitated  several 
jogs  of  up  to  3  km  from  the  profile  azimuth.  A  vehicle  odometer  was  used 
to  deploy  the  stations,  then  the  latitude  and  longitude  were  measured  using 
a  hand-held  GPS  receiver.  The  final  geographic  locations  listed  in  Table  2 
are  obtained  by  taking  multiple  RefTek  GPS  locations  at  each  site. 
Individual  GPS  locations  not  within  the  LI  -  sigma  of  the  median  location 
at  each  site  were  discarded.  The  remaining  GPS  locations  were  averaged 
to  give  a  final  location.  Final  geographic  locations  for  sites  that  either  did 
not  have  a  RefTek  GPS  installed,  or  did  not  receive  more  than  5  RefTek 
GPS  locations,  were  picked  off  a  USGS  30X60  minute  quadrangle  map. 

At  each  station,  the  sensors  were  placed  in  shallow  holes,  leveled  and  then 
covered  to  reduce  wind  noise.  The  recording  instruments  were  placed 
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several  feet  from  the  sensors  and  were  then  sheltered  from  the  sun  by  using 
garbage  bags  and  vegetation  cover  to  avoid  over-heating. 

3.3  Profile  Data 

The  data  from  both  Wyoming  Deep  Probe  explosions  have  been 
combined  to  form  a  seismic  profile  that  extends  from  about  134  km  to  280 
km  from  the  source  (Figure  3).  The  whole  profile  displays  exceptionally 
high  signal-to-noise  ratios  due  to  the  large  source  size  and  excellent 
coupling  in  the  quarry.  Unfortunately,  other  Deep  Probe  shots  outside 
Wyoming  were  of  only  marginal  quality  in  our  data  and  are  not  shown. 

The  complete  profile  data  set  is  currently  archived  in  SAC  format  on 
Exabyte  tape  and  is  currently  available  from  us  at  Phillips  Laboratory.  We 
intend  to  archive  the  data  at  the  IRIS  Data  Management  Center  at  a  later 
date. 


The  Wind  River  Mountains  -  Green  River  Basin  Seismic  Refraction 
Profile  can  be  used  for  a  variety  of  seismic  studies  due  to  similar  source 
properties  of  both  blasts.  We  verified  the  similarity  of  both  sources  by 
comparing  seismograms  from  each  shot  recorded  at  the  same  location. 
Stations  33  and  34  (Table  2)  were  occupied  for  the  first  and  second  blast, 
respectively.  Both  stations  are  at  the  same  location  and  both  were  equipped 
with  Refraction  Technology  recorders  and  HS-10  seismometers.  Figure  4a 
shows  time  domain  traces  from  the  two  blasts.  The  only  obvious 
difference  is  the  peak  at  108  seconds  in  blast  #1  (upper  panel)  which  is  not 
evident  in  blast  #2  (lower  panel).  Otherwise,  the  two  waveforms  are 
relatively  similar.  To  more  accurately  determine  the  differences  between 
the  two  blasts,  spectral  amplitudes  of  the  entire  blast,  the  P-wave,  and  a 
pre-event  ambient  noise  sample  are  examined  (Figures  4b  -  4d).  Figure  4b 
shows  the  spectral  amplitudes  of  a  six  minute  window  around  the  blasts. 

The  energy  of  both  blasts  is  predominantly  in  the  1.0  -  8.0  Hz  frequency 
band.  Within  this  frequency  band,  the  blasts  contain  different  spectral 
amplitude  peaks.  The  most  notable  is  that  blast  #1  (upper  panel)  contains  a 
higher  amplitude  peak  at  1.0  Hz  while  blast  #2  (lower  panel)  contains  a 
higher  peak  at  6.0  Hz.  Figure  4c  examines  the  spectral  amplitudes  of  only 
the  P-wave.  In  this  figure,  blast  #2  displays  a  more  prominent  peak  at  3.0 
and  7.0  Hz,  while  blast  #1  shows  a  peak  at  around  4.5  Hz.  Figure  4d 
displays  the  spectral  amplitudes  of  a  50  second  pre-event  noise  sample. 
Most  of  the  noise  occurs  between  0.0  and  1.0  Hz  (which  is  not  within  the 
1.0  -  8.0  Hz  frequency  band  in  which  the  seismic  signal  resides).  In 
addition,  the  seismic  noise  is  also  an  order  of  magnitude  smaller  in 
amplitude  than  the  seismic  signal.  The  preliminary  spectral  study  confirms 
that  although  the  seismic  signals  are  not  identical,  the  signals  from  both 
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blasts  have  enough  similarity  to  justify  combining  both  data  sets  for  this 
refraction  study. 

4.0  Future  Work 

This  report  documents  the  Wind  River  Mountains  -  Green  River 
Basin  Seismic  Refraction  Profile  data  set  using  2  large  chemical  explosions 
east  of  Lander,  Wyoming  as  the  source.  The  chemical  explosions  were 
recorded  across  the  Wind  River  Mountains  and  Green  River  Basin  of 
southwestern  Wyoming.  We  plan  to  use  this  data  set  to  constrain  the 
crustal  model  for  the  region  using  one  and  two  dimensional  travel  time  and 
waveform  modeling.  Finally,  the  seismograms  will  be  archived  for  use  by 
other  investigators. 
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Table  1. 


Shot  Locations 


<shnt#  Date  Time  (GMT) 

Latitude 

Longitude 

Elev.  (mV  Denth  fm)2 

143  08-09-95  11:30:00.000 

42.731  N 

107.667  W 

1958.0 

46.0 

243  08-17-95  11:30:00.000 

42.730  N 

107.665  W 

1958.0 

46.0 

Elevation  refers  to  quarry  water  level  relative  to  sea  level. 

2Depth  refers  to  distance  between  quarry  water  level  and  quarry  bottom  where  the 

ANFO  was  detonated. _ _ _ - — — - 
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Table  2.  Station  Locations 


Station  # 

Latitude  Longitude  Z  On) 

Recorder 

#  of  Channels 

Data  Available 

i 

42  32.09  N  109  12.54  W  2450 

TerraTek 

1 

No 

2 

42  33.84  N  109  14.85  W  2425 

TerraTek 

1 

No 

3 

42  34.20  N  109  17.14  W  2430 

TerraTek 

1 

Yes 

4 

42  33.13  N  109  19.45  W  2350 

TerraTek 

1 

Yes 

5 

42  33.58  N  109  21.80  W  2200 

TerraTek 

1 

No 

6 

42  36.73  N  109  24.11  W  2250 

PDAS 

3 

Yes 

7 

42  37.32  N  109  26.36  W  2250 

TerraTek 

1 

Yes 

8 

42  36.67  N  109  28.70  W  2190 

TerraTek 

1 

Yes 

9 

42  36.92  N  109  31.13  W  2180 

PDAS 

3 

Yes 

10 

42  35.42  N  109  33.50  W  2200 

PDAS 

3 

Yes 

11 

42  33.65  N  109  35.78  W  2215 

PDAS 

3 

Yes 

12 

42  36.22  N  109  38.20  W  2200 

PDAS 

3 

Yes 

13 

42  36.05  N  109  40.57  W  2210 

PDAS 

3 

Yes 

14 

42  36.17  N  109  42.88  W  2200 

PDAS 

3 

Yes 

15 

42  36.33  N  109  45.22  W  2150 

PDAS 

3 

Yes 

16 

42  35.95  N  109  48.64  W  2120 

PDAS 

3 

Yes 

17 

42  36.07  N  109  49.97  W  2100 

PDAS 

3 

Yes 

18 

42  35.99  N  109  52.18  W  2100 

PDAS 

3 

Yes 

19 

42  35.44  N  109  54.54  W  2100 

PDAS 

3 

Yes 

20 

42  33.95  N  109  56.90  W  2110 

PDAS 

3 

Yes 

21 

42  34.68  N  109  59.27  W  2110 

PDAS 

3 

Yes 

22 

42  34.67  N  110  01.53  W  2120 

PDAS 

3 

Yes 

23 

42  34.80  N  110  03.20  W  2110 

PDAS 

3 

Yes 

24 

42  34.07  N  110  06.24  W  2100 

PDAS 

3 

Yes 

25 

42  32.67  N  110  08.59  W  2100 

PDAS 

3 

Yes 

26 

42  32.59  N  110  10.96  W  2125 

PDAS 

3 

Yes 

27 

42  32.44  N  110  13.28  W  2150 

PDAS 

3 

Yes 

28 

42  32.37  N  110  15.62  W  2190 

RefTek 

3 

Yes 

29 

42  32.44  N  110  18.01  W  2250 

RefTek 

3 

Yes 

30 

42  31.92  N  110  20.35  W  2250 

RefTek 

3 

Yes 

31 

42  31.56  N  110  22.75  W  2300 

RefTek 

3 

Yes 

32 

42  31.61  N  110  24.99  W  2300 

RefTek 

3 

No 

33 

42  30.71  N  110  28.60  W  2380 

RefTek 

3 

Yes 

34 

42  30.71  N  110  28.60  W  2380 

RefTek 

3 

Yes 

35 

42  30.48  N  110  31.03  W  2450 

TerraTek 

1 

Yes 

36 

42  28.92  N  110  33.06  W  2500 

RefTek 

3 

Yes 

37 

42  27.50  N  110  35.29  W  2500 

TerraTek 

1 

No 

38 

42  28.82  N  110  37.50  W  2500 

RefTek 

3 

Yes 

39 

42  30.29  N  110  40.34  W  2600 

TerraTek 

1 

No 

40 

42  31.68  N  110  42.20  W  2650 

RefTek 

3 

No 

41 

42  31.63  N  110  44.65  W  2450 

TerraTek 

1 

Yes 

42 

42  29.69  N  110  48.12  W  2350 

RefTek 

3 

Yes 

43 

42  29.38  N  110  50.50  W  2300 

TerraTek 

1 

No 

44 

42  30.61  N  110  53.05  W  2250 

RefTek 

3 

Yes 

45 

42  29.17  N  110  55.12  W  2150 

TerraTek 

1 

Yes 

46 

42  23.99  N  111  00.46  W  1950 

RefTek 

3 

Yes 

47 

42  24. 14  N  11 1  02.70  W  1950 

TerraTek 

1 

Yes 
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Table  2.  (Continued) 


Mema  Jet  42  56.34  N  1 10  20.79  W  2340 
Big  Sandy  42  37.91  N  109  28.04  W  2200 
Big  Piney  42  32.06  N  110  16.53  W  2250 
Fontenelle  42  05.44  N  1 10  10.06  W  2000 
Allred  42  29.55  N  110  57.74  W  1900 


RelTek 

RefTek 

RefTek 

RefTek 

RefTek 


3 

3 

3 

3 

3 


Yes 

Yes 

Yes 

Yes 

Yes 
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List  of  Figures 


FIGURE  1. 


FIGURE  2. 


FIGURE  3. 


FIGURE  4a. 


FIGURE  4b. 


FIGURE  4c. 


FIGURE  4d. 


Location  map.  The  Green  River  refraction  profile  is  shown 
by  the  heavy  dashed  line.  The  Wyoming  shot  points  are 
denoted  by  the  star  SE  of  Riverton.  The  profile  was  selected 
to  take  advantage  of  large-aperture  array  stations  (house 
symbol).  Dotted  lines  indicate  approximate  location  of  the 
Wind  River  and  Wyoming  Range  Mountains. 


Refraction  Station  Schematic.  Stations  consisted  of  sensor, 
recorder,  battery  power,  and  timing  system.  The 
arrangement  for  a  Refraction  Technology  recorder  station  is 
shown  here. 


Green  River  Seismic  Profile  record  section.  Data  from  both 
shots  are  displayed  in  trace  normalized  record  section 
format  reduced  at  6  km/sec.  East  (Big  Sandy)  is  to  the  left; 
the  right  side  is  at  the  Idaho-Wyoming  border. 


Seismograms  from  co-located  stations  for  shot  143  (upper 
panel)  and  243  (lower  panel).  Windows  for  noise  and  P- 
wave  spectra  are  shown  at  bottom. 


Amplitude  spectra  for  entire  window  shown  in  Figure  4a. 
Most  of  the  energy  is  between  0.5  and  8  Hz,  although 
spectral  peaks  are  different. 


Amplitude  spectra  for  P-wave  window  (see  Figure  4a).  Main 
peaks  correlate  in  frequency,  but  not  in  amplitude.  Note 
extra  energy  at  ~5  Hz  for  shot  143. 


Amplitude  spectra  for  noise  window  (see  Figure  4a).  Noise  is 
predominantly  below  1  Hz  and  is  a  factor  of  200  less  than 
the  signal  in  the  2  -  4  Hz  range. 
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